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Abstract 

This  thesis  discusses  the  recent  developments  to  the  electric  propulsion  plume 
code  DRACO  as  well  as  a  validation  and  sensitivity  analysis  of  the  code  using  data  from 
an  AFRL  experiment  using  a  Busek  200  W  Hall  Thruster.  DRACO  is  a  PIC  code  that 
models  particles  kinematically  while  using  finite  differences  schemes  to  solve  the  electric 
potential  and  field. 

The  DRACO  code  has  been  recently  modified  to  improve  simulation  results, 
functionality  and  performance.  A  particle  source  has  been  added  that  uses  the  Hall 
Thruster  device  code  HPHall  as  input  for  a  source  to  model  Hall  Thrusters.  The  code  is 
now  also  capable  of  using  a  non-uniform  mesh  that  uses  any  combination  of  uniform, 
linear  and  exponential  stretching  schemes  in  any  of  the  three  directions.  A  stretched  mesh 
can  be  used  to  refine  simulation  results  in  certain  areas,  such  as  the  exit  of  a  thruster,  or 
improve  performance  by  reducing  the  number  of  cells  in  a  mesh.  Finally,  DRACO  now 
has  the  capability  of  using  a  DSMC  collision  scheme  as  well  as  performing 
recombination  collisions. 

A  sensitivity  analysis  of  the  newly  upgraded  DRACO  code  was  performed  to  test 
the  new  functionalities  of  the  code  as  well  as  validate  the  code  using  experimental  data 
gathered  at  AFRL  using  a  Busek  200  W  Hall  Thruster.  A  simulation  was  created  that 
attempts  to  numerically  recreate  the  AFRL  experiment  and  the  validation  is  performed  by 
comparing  the  plasma  potential,  polytropic  temperature,  ion  number  density  of  the 
thruster  plume  as  well  as  Faraday  and  ExB  probe  results.  The  study  compares  the  newly 
developed  HPHall  source  with  older  source  models  and  also  compares  the  variations  of 
the  HPHall  source.  The  field  solver  and  collision  model  used  are  also  compared  to 
determine  how  to  achieve  the  best  results  using  the  DRACO  code.  Finally,  both  uniform 
and  non-uniform  meshes  are  tested  to  determine  if  a  non-uniform  mesh  can  be  properly 
implemented  to  improve  simulation  results  and  performance. 

The  results  from  the  validation  and  sensitivity  study  show  that  the  DRACO  code 
can  be  used  to  recreate  a  vacuum  chamber  simulation  using  a  Hall  Thruster.  The  best 
results  occur  when  the  newly  developed  HPHall  source  is  used  with  a  MCC  collision 
scheme  using  a  projected  background  neutral  density  and  CEX  collision  tracking.  A 
stretched  mesh  was  tested  and  proved  results  that  are  as  accurate  as  a  uniform  mesh,  if 
not  more  accurate  in  locations  of  high  mesh  refinement. 
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Chapter  1 :  Introduction 


The  research  of  this  thesis  focuses  on  using  the  plasma  simulation  code 
COLISEUM  to  observe  the  plume  of  a  Hall  Thruster  in  a  vacuum  chamber.16  To  validate 
the  code,  the  results  of  the  simulations  are  compared  to  experimental  results.  In  this 
chapter  the  basics  of  electric  propulsion  are  outlined  including  the  interactions  between 
the  spacecraft  and  the  thruster  plasma  environment.  The  COLISEUM  code  is  also  briefly 
introduced.  Finally,  the  remainder  of  the  thesis  is  outlined. 

1.1 :  Basics  of  Electric  Propulsion 

Electric  propulsion  is  an  innovative  method  of  producing  low  amounts  of  thrust  to 
perform  large  changes  in  velocity  over  a  long  period  of  time.  The  thrusts  produced  are 
too  small  to  send  a  spacecraft  into  orbit;  however,  once  in  orbit  a  spacecraft  can  use 
electric  propulsion  systems  to  increase  or  maintain  its  orbit  radius.  The  changes  in 
velocity  can  be  orders  of  magnitude  greater  than  that  of  conventional  chemical  propulsion 
systems  however  the  low  thrust  often  requires  a  longer  travel  time  for  the  same  propellant 
mass. 


There  are  several  different  electric  propulsion  systems  currently  available.  Two  of 
the  most  common  types  are  Ion  Thrusters  and  Hall  Thrusters.  These  two  propulsion 
systems  are  similar  in  that  they  both  create  plasma  and  then  use  electrostatic  forces  to 
accelerate  the  positively  charged  ions  in  the  plasma  to  produce  thrust.  In  Ion  thrusters  the 
ions  are  accelerated  by  ion  optics,  a  series  of  charged  flat  plates  with  small  holes  for  the 
ions  to  travel  through.  In  Hall  thrusters  the  ions  are  accelerated  by  the  JxB  force.  In  both 
cases,  after  the  ions  are  accelerated,  they  are  then  neutralized  with  an  electron  beam  to 
prevent  the  ions  from  returning  to  the  spacecraft.  The  research  in  this  thesis  focuses  on 
Hall  Thrusters. 

A  major  component  in  the  plume  of  an  Ion  Thruster  or  Hall  Thruster  is  the 
charge-exchange  ions  generated  by  CEX  collisions  between  beam  ions  and  neutrals.  It  is 
well  known  that  CEX  ions  can  be  pushed  out  of  the  plume  and  back  flow  to  dominate  the 
plasma  environment  surrounding  the  spacecraft.28  The  backflow  of  CEX  ions  can  also 
cause  ion  sputtering  from  thruster  components.  These  sputtered  ions  from  thruster 
components  can  cause  spacecraft  contamination.8  The  interactions  between  the  spacecraft 
and  the  plasma  environment  created  by  the  electric  thrusters  must  be  accurately  studied  to 
avoid  mission  failures. 

1.1.1:  Hall  Thrusters 

The  governing  principal  behind  a  Hall  Thruster  is  the  Hall  Effect,  which  traps 
electrons  in  a  cloud  at  the  exit  of  the  thruster.  The  Hall  Effect  was  discovered  in  1879  by 
the  American  physicist  Edwin  Hall.13  The  Hall  Effect  is  a  potential  gradient  caused  by 
charged  particles  traveling  normal  to  an  electric  field  which  is  also  perpendicular  to  a 
magnetic  field.  The  magnetic  field  is  produced  by  two  electro-magnets,  one  outer 
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magnetic  coil  around  the  outer  ring  of  the  discharge  channel  and  one  inner  magnet  around 
the  central  axis  of  the  thruster.  Electrons  are  emitted  from  an  external  cathode  and  are 
attracted  into  the  thruster  channel  and  the  particles  then  revolve  around  the  central  axis  of 
the  thruster.  The  electrons  remain  trapped  by  the  Hall  Effect  at  the  exit  of  the  thruster 
where  the  magnetic  field  is  at  its  greatest  strength.  The  electrons  trapped  in  this  cloud 
have  a  high  energy  (10-20eV).  When  neutral  gas  is  injected  from  the  thruster  channel  the 
neutral  particles  are  bombarded  by  the  high  energy  electrons  and  become  ionized.  The 
ions  are  accelerated  by  the  electric  potential  created  between  the  anode  at  the  back  of  the 
thruster  channel  and  the  external  cathode,  producing  thrust.13  A  schematic  of  a  typical 
Hall  Thruster  is  shown  below  in  Figure  1.1. 


Boron  Nitride  Walls 


Anode/Gas 

Distributor 


Cathode  Neutralizer 


Inner 

Magnetic 

Coil 


Magnetic  Circuit  Outer  Magnetic  Coil 


Figure  1.1:  Section  Schematic  of  Hall  Effect  Thruster11 


The  propellant  most  commonly  used  for  Hall  Thrusters  is  Xenon  because  of  its 
high  molecular  weight  and  low  ionization  potential.  Most  of  the  Xenon  particles  are 
singly  charged  but  about  10%  of  the  particles  are  doubly  charged.  Additionally,  the 
ionization  process  is  about  90%  efficient  so  some  neutral  gas  is  also  being  emitted  along 
with  the  ions.12  CEX  collisions  will  occur  between  the  fast  moving  beam  ions  and  slow 
moving  neutral  atoms  creating  slow  moving  CEX  ions  and  fast  moving  neutrals.  The 
typical  discharge  voltage  is  approximately  200-500V  and  the  typical  specific  impulse  of  a 
Hall  Thruster  is  approximately  1000-2000s.9 

The  thruster  that  is  the  focus  of  this  thesis  is  a  Busek  200W  Hall  Thruster.  A 
photograph  of  the  thruster  is  shown  in  Figure  1.2  and  this  thruster  has  been  extensively 
tested  at  the  Air  Force  Research  Laboratory  (AFRL)  at  Edwards  Air  Force  Base  in 
California.  Figure  1.3  is  a  photograph  of  the  200W  Hall  thruster  firing  in  a  vacuum 
chamber  at  AFRL. 
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Figure  1.3:  Photograph  of  Busek  200 W  Hall  Thruster  Firing  in  Vacuum  Chamber  at  Air  Force 

Research  Laboratory 

1.1.2:  Spacecraft  -  Environment  Interactions 

The  plume  environment  created  by  an  electric  propulsion  (EP)  device  has  long 
been  a  major  concern.  The  effects  of  EP  plume  range  from  plume  contamination,  plasma 
interaction  effects  such  as  charging  and  interference  with  solar  array,  to  interference  with 
plasma  measurements.  Obviously,  these  effects  have  a  major  impact  on  a  mission  using 
an  electric  propulsion  system  and  must  be  studied  for  each  mission. 

1.2:  Modeling  Electric  Propulsion  Plume  Interactions 

Numerical  modeling  has  been  increasingly  used  to  study  the  complex  plasma 
behavior  in  an  EP  thruster.  For  instance,  PIC  code  developed  by  Wang  et.  al.27  was  used 
to  predict  the  Ion  Thruster  plume  for  the  Deep  Space  1  Spacecraft  and  the  results  were 
found  to  be  in  excellent  agreement  with  In-Flight  measurements  from  Deep  Space  1 .  This 
thesis  concerns  with  Hall  Thruster  plume.  Some  recent  work  concerning  Hall  Thruster 
plume  modeling  is  summarized  in  Boyd  5. 
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1.3:  COLISEUM 


The  numerical  code  used  in  this  thesis  is  COLISEUM,  a  plasma  simulation 
framework  developed  at  the  Air  Force  Research  Laboratory  at  Edwards  Air  Force  Base  in 
collaboration  with  Virginia  Tech  and  the  Massachusetts  Institute  of  Technology.16  The 
COLISEUM  framework  handles  input  and  output  for  all  of  the  COLISEUM  modules. 
The  COLISEUM  simulation  modules  are  DRACO,7  developed  at  Virginia  Tech, 
AQUILA,  developed  at  MIT,  and  RAY.  Each  module  solves  a  different  type  of  plasma 
simulation. 

DRACO  uses  the  Particle-in-Cell  (PIC)  algorithm2  to  solve  plasma  simulation 
problems.  PIC  is  a  method  for  numerically  solving  plasma  simulations  that  uses  a  kinetic 
description  of  plasma  with  the  electric  field  solved  on  a  mesh.  The  plasma  is  tracked  on  a 
Cartesian  mesh  and  surface  boundaries  are  defined  using  a  triangular  mesh.  Particle 
collisions  can  be  handled  using  the  Monte  Carlo  and  Direct  Simulation  Monte  Carlo 
methods.  DRACO  is  designed  to  work  with  several  different  field  solvers  including  both 
Quasi-Neutral  and  Non-Neutral  field  solvers.  DRACO  is  the  module  used  in  all  the 
research  for  this  thesis.  A  more  detailed  discussion  of  the  DRACO  module  is  provided  in 
the  next  chapter. 

The  other  two  COLISEUM  modules  are  AQUILA  and  RAY.  AQUILA  is  a 
hybrid-PIC  code  that  uses  a  body-fitted,  unstructured  tetrahedral  mesh.  It  also  handles 
particle  collisions  using  the  Direct  Simulation  Monte  Carlo  method.24  RAY  is  a 
sputtering  analysis  module  that  uses  ray  tracing  to  determine  where  sputtering  will  occur 
and  where  sputtered  material  will  be  deposited  on  a  spacecraft.15 

1.4:  Thesis  Overview 

The  primary  objective  of  this  thesis  is  to  use  the  DRACO  module  of  COLISEUM 
to  simulate  the  plasma  plume  of  a  Hall  Thruster  in  a  vacuum  chamber.  The  purpose  of 
this  research  is  to  validate  the  DRACO  code  by  comparing  solution  results  to 
experiments  conducted  at  AFRL.  By  using  DRACO,  any  user  can  gain  knowledge  about 
how  a  thruster  plume  will  act  before  spending  money  on  vacuum  chamber  tests. 
Additionally,  DRACO  can  be  used  to  conduct  real  in-space  simulations  and  create 
environments  that  are  impossible  for  experiments. 

Chapter  2  of  this  thesis  will  provide  a  more  detailed  description  of  COLISEUM, 
in  particular  the  DRACO  module  which  is  the  focus  of  this  research.  The  different 
algorithms  used  by  DRACO  including  the  equations  used  by  the  program  are  outlined  in 
detail.  In  Chapter  3  the  newly  developed  functions  of  DRACO  are  described.  This 
includes  a  new  source  model  using  the  Hall  Thruster  device  code  HPHall,  DSMC 
capabilities  for  DRACO  including  functionality  for  recombination  collisions  and  the 
ability  to  use  a  non-uniform  Cartesian  mesh.  In  Chapter  4  the  set-up  for  the  simulations  is 
described  as  well  as  the  experiments  conducted  at  AFRL.  In  Chapter  5  the  results  of  all  of 
the  test  cases  are  given  including  plots  and  tables.  In  Chapter  6,  conclusions  are  drawn  as 
well  as  suggestions  for  future  research  work. 
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Chapter  2:  Numerical  Methods 


DRACO  is  a  C  code  originally  developed  from  Dr.  Joseph  Wang’s  Fortran  77 
PLUME  code.27  The  code  is  part  of  the  COLISEUM  infrastructure  and  is  capable  of  both 
Full-PIC  and  Hybrid-PIC  simulations.  DRACO  is  equipped  with  multiple  potential  field 
solvers  as  well  as  a  series  of  particle  collision  models.  The  meshing  is  handled  by  the 
DRACO  module  YOLCAR.  The  major  difference  between  DRACO  and  the  other 
modules  is  that  DRACO  used  a  Cartesian  based  PIC  mesh  and  supports  two  formulations, 
the  standard  finite-difference  and  the  newly  developed  Immersed  Finite  Element 
formulation. 

2.1 :  Governing  Equations 

A  charged  particle  experiences  a  force  in  the  presence  of  an  electric  field  as  well 
as  a  force  equal  to  the  cross  product  of  the  velocity  of  the  particle  and  a  magnetic  field. 
This  is  known  as  the  Lorentz  force  and  is  the  governing  force  on  all  particles  in  a  plasma 
plume  simulation.29 

F  =  q[E  +  vxB)  (2.1) 

The  Lorentz  force,  F ,  is  a  function  of  the  particle  charge,  q,  electric  field,  E , 

particle  velocity,  v  ,  and  magnetic  field,  B  .  The  COLISEUM  code  is  an  electrostatic  PIC 
code  and  assumes  that  the  plasma  induced  magnetic  field  is  negligible.  Thus,  the 
magnetic  field  is  a  defined  function  and  usually  assumed  to  be  constant.  If  the  electric 
and  magnetic  field  are  known  at  any  given  time  then  the  forces  on  the  particle  can  be 
determined  assuming  that  the  Lorentz  force  is  the  only  force  on  the  particle.6  In  addition 
to  the  Lorentz  forces  on  the  particles,  particle  and  surface  collisions  can  also  be  tracked  in 
order  to  fully  describe  the  motion  of  the  particles.  The  electric  field  must  be  solved  to 
determine  the  Lorentz  force. 


V-E  =  ^-  (2.2) 

fco 

The  gradient  of  the  electric  field  is  equal  to  the  charge  density,  p,  divided  by  the 
permittivity  of  free  space,  e0 .  The  electric  field  can  also  be  related  to  the  scalar  electric 
potential,  </> . 


E  =  -V</>  (2-3) 

The  scalar  electric  potential  can  be  solved  using  the  Poisson’s  equation  and  then 
the  electric  field  can  be  solved  using  the  potential. 

VV  =  -t-  (2.4) 

t0 
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2.2:  Finite  Difference  Method 


Finite  difference  methods  are  used  to  solve  the  governing  equations  for  the 
electric  potential  and  electric  field.  A  central  finite  difference  scheme  is  used  for  solving 
the  governing  equations.  The  central  finite  difference  scheme  is  used  for  approximating 
the  derivatives.  The  domain  of  a  simulation  is  divided  into  a  mesh  where  the  potential 
and  electric  field  are  calculated  at  each  grid  point  in  the  mesh.  The  central  difference 
approximation  for  the  first  derivative  of  a  function  is  given  in  Equation  2.5.  Similarly  the 
central  difference  approximation  for  the  second  derivative  of  a  fuction  is  given  in 
Equation  2.6.  These  equations  are  calculated  at  a  node  and  use  the  value  at  the  node, 
fi  j  k ,  as  well  as  the  value  of  the  node  ahead,  fi+ljk ,  and  behind,  /(_l  ;  /{ ,  that  node.  These 

equations  also  use  the  distance  between  the  current  node  and  previous  node,  Av_ ,  and  the 
distance  between  the  current  node  and  the  next  node,  Ax+ . 


3/  +(aV  -  a*;  +  (-A r:  t 

dx  AxlAx_  +  Ax2_Ax+ 

32/  _  2(Ai-+  )f,_]jk  +  2(- Ax_  -Av+)/.m  +2(-A x_)fMJ<k 
dx2  Ax^Ax_  +  Ax2Ax+ 


If  a  uniform  mesh  is  used  then  these  equations  can  be  simplified  because  Ax_  is 
equal  to  Av+ .  The  simplified  versions  of  these  equations  are  shown  below. 


df  _  f i i+l,j,k  fi-\,j,k 

dx  2Ax 


(2.7) 


d2/  _  f i i+l,j,k  ^fi,j,k  + 

~d^~  Ax2 


(2.8) 


The  central  difference  approximation  for  a  second  derivative  is  plugged  into  the 
Poisson’s  equation  to  solve  for  the  potential.  The  approximation  is  made  for  all  three 
directions  and  modifications  to  the  code  have  been  made  to  accommodate  for  non- 
uniform  meshes  as  well.  Further  discussion  of  these  code  developments  are  discussed  in 
detail  later  in  this  thesis. 


2.3:  Particle  Motion 


In  DRACO,  the  trajectories  of  each  particle  are  solved  from  Newton’s  2nd  Law. 
The  position  and  velocity  of  the  particles  are  calculated  at  a  given  time  using  the  previous 
time  steps  position  and  velocity  and  Equation  2.1  to  determine  the  force  on  the  particle. 
The  particle’s  charge  is  deposited  onto  the  Cartesian  mesh  based  on  its  position  in  a  cell. 
A  detailed  description  of  the  particle  motion  in  DRACO  is  discussed  in  Reference  5. 
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2.4:  Particle  Collisions 


An  important  aspect  of  a  Particle-in-Cell  code  is  tracking  collisions  between 
particles.  The  two  typical  methods  for  handling  particle  collisions  in  PIC  codes  are  Monte 
Carlo  Collisions  (MCC)  and  Direct  Simulation  Monte  Carlo  (DSMC).  DRACO  has  been 
capable  of  performing  MCC  collisions  for  some  time  and  DSMC  functionally  has  been 
recently  added  and  will  be  discussed  in  Chapter  3.  The  collisions  that  DRACO  is  capable 
of  are  Charge  Exchange  (CEX),  Variable  Hard  Sphere  (VHS)  and  recombination 
collisions. 

2.4.1:  Monte  Carlo  Collisions 

A  Monte  Carlo  Collision  is  numerically  performed  by  treating  the  ion  particles  as 
simulated  particles  and  the  neutral  particles  as  a  background  neutral  gas.  This  method 
reduces  the  computational  demands  of  the  particle  collisions  because  particles  do  not 
need  to  be  individually  paired  with  each  other  to  determine  if  a  collision  has  occurred. 
There  are  two  major  steps  in  performing  a  particle  collision:  determine  if  a  particle 
collision  has  occurred  and  then  perform  the  corresponding  collision.  A  collision  cross- 
section,  <j ,  is  used  to  determine  if  a  particle  has  collided  with  another  particle. 

2.4.2:  Charge  Exchange 

Charge  exchange  collisions  occur  when  a  fast  moving  charged  particle  collides 
with  a  slow  moving  neutral  particle.  The  charge  from  the  ion  is  transferred  into  the 
neutral  particle  and  the  result  is  a  fast  moving  neutral  particle  and  a  slow  moving  ion.  The 
collision  cross-section  is  dependant  on  the  relative  velocity  between  particles,  cr .  The 
velocity  of  the  neutral  particles  is  assumed  to  be  zero  in  this  analysis.  There  are  multiple 
collision  cross-section  models  for  each  type  of  particle  collision.  For  this  thesis  the  cross- 
section  used  was  developed  by  Pullins  and  scales  with  the  log  of  the  relative  velocity  of 
the  particles.25 


<JCEXP  =1.1872xlO"20[-23.31og(cr)  + 188.81]  (2.9) 

Using  the  collision  cross-section,  the  probability  of  a  collision,  P,  is  computed 
using  the  particle  number  density,  n,  and  time  step,  At.  The  probability  equation  produces 
a  number  between  zero  and  one  and  a  random  number  is  chosen  between  zero  and  one  to 
determine  if  a  particle  collision  occurs.  If  the  random  number  is  greater  than  the 
probability  then  a  particle  collision  is  performed.25 

P  =  [l.0-exp(-ncr  crAf)]  (2.10) 

The  resulting  ion  produced  by  a  CEX  collision  has  a  velocity  based  on  the 
Maxwellian  distribution  of  the  thermal  velocity  of  the  neutral  gas  in  a  random  direction. 
The  thermal  velocity  is  based  on  the  Boltzmann  constant,  k,  gas  temperature,  T,  and 
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particle  mass,  m.  The  Maxwellian  distribution  is  based  on  M  random  numbers,  and  for 
the  CEX  collisions  performed  by  DRACO,  3  random  numbers  are  used. 


'  [2kT^ 

M 

Af" 

~M~ 

m  J 

HRi 

_  i= 1 

T_ 

12 

2.4.3:  Variable  Hard  Sphere 


(2.11) 


Variable  Hard  Sphere  (VHS)  collisions  are  elastic  collisions  between  two 
particles.  VHS  collisions  must  be  handled  differently  than  CEX  collisions  because 
particles  must  be  paired  together  to  conserve  momentum.  To  do  this  a  hybrid  MCC 
method  is  used  to  group  the  particles  into  cells  and  then  pair  each  particle  with  a  random 
target  particle.  The  probability  of  a  VHS  collision  is  orders  of  magnitude  lower  than  that 
of  CEX  collisions.  The  collision  cross-section  used  for  calculating  the  collision 
probability  for  VHS  collision  is  shown  below. 


®VHS 


8.2807x10 


-16 


(2.12) 


The  probability  of  collision  is  calculated  using  the  VHS  collision  cross-section 
and  the  relative  velocity  between  the  source  and  target  particles.  The  probability  is 
calculated  in  the  same  way  as  a  CEX  collision  and  compared  to  a  random  number.  To 
perform  the  collision,  random  angles,  %  and  e ,  are  chosen  between  the  two  particles  and 
an  elastic  collision  is  performed  along  that  angle.  The  angle  X  is  defined  below  and  the 
angle  e  is  a  random  number  between  zero  and  2 n  ,25 

X  -  cos~l(2rnd  -l)  (2.13) 


The  next  step  is  to  calculate  the  center  of  mass  velocity,  cm ,  using  the  masses  and 
velocities  of  the  source  and  target  particles.  Using  the  random  angles  and  the  center  of 
mass  velocity  the  collision  adjustment  velocity,  (O ,  can  be  calculated.25 


mc  m, 

Cm= - - — vs+ - - — v,  (2.14) 

ms  +mt  ms  +mt 

a>  =  cr  cos^)*  +  cr  sin(/jf)cos(£  )j  +  c  r  sin (x) sin (e)k  (2.15) 

The  final  velocities  of  the  collided  source  and  target  particles  can  be  determined 
using  the  particle  collision  adjusted  velocity,  particle  masses  and  the  center  of  mass 
velocity.25 
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G 

~Cm  + 

CO 

(2.16) 

ms  +  mt 

m, 

G 

~Cm  + 

CO 

(2.17) 

ms  +  mt 

2.5:  Field  Solver 


The  correct  way  to  solve  the  potential  in  an  electrostatic  PIC  simulation  is  to  solve 
of  the  Poisson’s  equation  on  a  Cartesian  mesh.  The  Poisson’s  equation  can  be  solved  by 
more  than  one  method.  A  fully-kinetic  plasma  simulation  can  be  run  where  both  ions  and 
electrons  are  tracked  but  this  causes  the  simulation  to  be  very  computational  demanding. 
A  quasi-neutral  simulation  can  be  run  by  assuming  that  the  electron  density  is 
approximately  the  same  as  the  ion  density.  This  removes  the  need  to  track  electrons  but 
requires  that  a  non-linear  solver  must  be  used  to  solve  the  Poisson’s  equation.  DRACO  is 
equipped  with  several  non-linear  solvers:  Gauss-Seidel,  Dynamic  Alternate  Direct 
Implicit  (DADI)  and  Preconjugate  Gradient  (PCG). 


If  certain  assumptions  are  made  then  the  potential  can  be  calculated  by  means  of 
the  Boltzmann  inversion  of  the  electric  field  equation.  The  assumptions  are  that  the 
simulations  have  high  plasma  densities  and  the  plasma  is  quasi-neutral.  Additionally,  the 
Boltzmann  solver  can  not  resolve  the  plasma  sheath.10  This  method  does  not  use  the 
Poisson’s  equation  but  instead  assumes  that  the  potential  is  directly  related  to  charge 
density.  This  method  greatly  simplifies  the  field  solving  and  is  a  very  fast  and  simple 
approach.  The  emphasis  of  this  research  is  to  compare  results  in  the  center  plume  region 
where  the  plasma  is  quasi-neutral.  Hence,  to  save  computational  time,  most  of  the 
simulations  in  this  thesis  used  a  simplified  approach:  we  assume  that  the  electrons  can  be 


approximated  by  the  Boltzmann  distribution,  ne  —  n0  exp 


and  we  solve  for  the 


V  o  j 

potential  directly  from  rii=ne.  Note  this  approach  is  only  valid  for  the  center  plume  region 
and  cannot  resolve  the  plasma  sheath. 


2.5.1 :  Boltzmann  Solver 


The  Boltzmann  solver  is  able  to  calculate  the  potential  field  using  the  charge 
density,  p ,  calculated  during  the  particle  push  and  either  a  constant  or  polytropic 

temperature  model.  The  solver  requires  a  user  defined  reference  potential,  (j){] , 
temperature,  T0 ,  number  density,  n0 ,  and  specific  heat  ratio,  y .  The  equation  for 
potential  for  a  given  node,  n,  is  listed  below  for  a  constant  temperature  model. 


<t>n  ~  ^0  + 


kTn 


-In 


(  p  ^ 

\eno  J 


(2.18) 
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A  polytropic  temperature  model  can  be  used  with  the  Boltzmann  inversion  to 
calculate  the  potential  as  well. 


r  ) 

- 1 

f  Pn  1 

Y- 1 

kT0 

lr-ij 

e 

Ken0  J 

e 

(2.19) 


The  Boltzmann  solver  is  the  field  solver  used  in  this  thesis  because  the 
simulations  conducted  are  all  near  the  exit  of  the  Hall  Thruster  and  have  high  densities  so 
the  assumptions  of  the  Boltzmann  inversion  still  apply. 
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Chapter  3:  Code  Development 


The  DRACO  code  has  undergone  several  major  upgrades  in  the  past  year 
including  the  addition  of  non-uniform  mesh  capability,  Direct  Simulation  Monte  Carlo 
(DSMC)  collision  handling,  capability  for  recombination  collisions  and  the  addition  of  a 
new  source  model  based  on  the  device  code  HPHall.  The  DSMC  code  was  added  by 
incorporating  the  Collide  library  developed  by  Professor  Boyd  at  the  University  of 
Michigan.  The  recombination  code  was  written  using  the  Collide  library  and  the  new 
HPHall  source  was  created  using  a  previous  DRACO  source  as  a  template.  The  basic 
functionality  and  process  of  the  code  development  is  discussed  below.  The  non-uniform 
code  development,  integration  of  the  DSMC  code  and  writing  of  the  recombination  code 
was  done  in  collaboration  with  Lubos  Brieda  and  Alex  Barrie. 

3.1 :  Non-Uniform  Mesh 

Simulations  involving  a  large  domain  and  a  large  number  of  particles  have  a  long 
processing  time,  often  days.  A  domain  which  implements  a  stretched  mesh,  varying  cell 
size,  can  be  used  to  reduce  memory  requirements  because  fewer  cells  are  required.  The 
DRACO  code  was  recently  developed  to  use  a  non-uniform  mesh  with  both  linear  and 
exponential  stretching  schemes.  A  stretched  mesh  has  the  advantage  of  having  a  fine  cell 
mesh  near  the  exit  of  a  thruster  where  the  densities  are  high  and  a  coarse  mesh  further 
downstream  of  the  plume  where  densities  are  lower.  The  stretched  meshes  are  defined  in 
the  domain  input  file  with  an  initial  cell  size,  Axo,  the  number  of  cells  in  the  stretched 
zone,  n,  initial  position,  xo,  and  final  position  of  the  zone,  xmax.  The  code  calculates  a 
stretching  factor,  k,  and  the  node  positions  based  on  the  type  of  stretched  mesh  used.  The 
one  dimensional  equations  for  the  position  of  the  nodes  for  both  the  linear  (Equation  3.1) 
and  exponential  (Equation  3.3)  stretching  schemes  are  listed  below.  Additionally  the  cell 
size,  dx,  must  be  determined  for  each  cell  (Equation  3.2,  3.4)  because  it  will  vary  for  each 
cell  and  is  needed  for  the  potential  and  electric  field  equations. 

Linear: 


Exponential: 


x  — 


A x0i 

Vl-0.5£y 


+ 


Ax0k 

vl-0.5  kj 


1  -2  • 
—  I  -l 
2 


+  xn 


dx  ■ 


Axn 


(1-0.5  k) 


(l  +  ki) 


e  -1 

dx  =  -^-(ek{M)  -eu) 
ek -1 


(3.1) 

(3.2) 


(3.3) 

(3.4) 
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The  next  step  of  implementing  a  stretched  mesh  is  being  able  to  determine  the 
node  to  which  the  particle  belongs.  This  is  preformed  using  Equation  3.5  for  a  linear 
mesh  and  Equation  3.6  for  an  exponential  stretched  mesh  and  rounding  the  values  down 
to  the  nearest  integer.  This  functionality  is  used  by  DRACO  during  scatter  and  gather 
operations  to  determine  where  to  deposit  the  charge  of  a  particle  onto  the  mesh. 


Linear: 


Exponential: 


i  = 


(k-l)+k 


Axn 


W 


l -0.5k 


(1  ~k) 


-2 


{x0  -  x ) 


(3.5) 


(3.6) 


The  final  modification  to  the  code  that  is  required  in  order  to  implement  a 
stretched  mesh  is  to  modify  the  scatter  and  gather  operations  as  well  as  the  finite 
difference  equations  to  incorporate  a  stretched  mesh.  Originally  the  code  was  developed 
assuming  that  the  cell  size  is  uniform  so  all  cell  size  references,  Ax ,  must  be  replaced 
with  the  corresponding  left  and  right  cell  size  for  each  node,  Ax_  and  Ax+ .  There  are 
several  instances  for  which  the  code  must  be  modified  with  this  change  and  the  code  was 
updated  accordingly.  After  the  code  had  been  modified  a  test  case  was  run  using  three 
domains  of  equal  size  with  three  different  mesh  types:  uniform,  linear  stretch  and 
exponential  stretch.  The  steady  state  potential  solution  for  the  flow  around  a  sphere  is 
shown  below  for  all  three  meshes.  All  three  meshes  have  a  uniform  mesh  in  the  x  and  y 
directions  with  the  same  number  of  cells  for  each  case  but  the  uniform  mesh  has  30  nodes 
in  the  z-direction  and  the  non-uniform  meshes  have  20  nodes  in  the  z-direction.  A  2D 
contour  plot  of  a  slice  through  the  center  of  the  sphere  for  all  three  meshes  is  shown  in 
Figure  3.1  -  Figure  3.3. 

It  can  be  seen  in  these  plots  that  the  results  of  the  stretched  mesh  are  almost 
exactly  the  same  for  all  three  cases.  Some  resolution  is  lost  in  the  largest  cells  but  the 
results  in  the  fire  grid  region  are  identical  for  all  three  cases.  Although  some  resolution  is 
lost  in  the  large  cells,  this  occurs  in  low  density  regions  in  the  wake  of  the  sphere  where 
high  resolution  is  not  required.  Further  analysis  of  the  non-uniform  mesh  capability  of 
DRACO  is  preformed  as  part  of  this  thesis  and  is  discussed  in  detail  in  Chapter  5. 
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3.2:  Direct  Simulation  Monte  Carlo 


Direct  Simulation  Monte  Carlo  (DSMC)  is  a  method  for  tracking  particle 
collisions  in  a  computational  simulation  that  collides  two  actual  macro  particles,  unlike 
MCC  which  collides  one  macro  particle  with  a  background  gas.1  Particles  are  grouped 
into  collision  pairs  and  the  probability  of  collision  for  each  pair  is  computed.  A  DSMC 
simulation  tends  to  be  more  computationally  demanding  because  collisions  involving 
neutral  particles  must  track  the  neutral  particles  rather  than  simply  pre-computing  the 
background  neutral  densities.  In  addition  to  greatly  increasing  the  number  of  particles  in 
the  simulation  the  code  must  also  track  what  cells  the  particles  reside  in  at  any  given  time 
so  that  the  particles  can  be  grouped  into  collision  pairs.  Simulations  using  DSMC  rather 
than  MCC  tend  to  take  longer  but  produce  results  that  are  more  physical  because  DSMC 
insures  that  momentum  is  conserved  in  the  simulation. 

DRACO  was  recently  upgraded  for  DSMC  capability  by  use  of  the  collide  library 
developed  by  Iain  Boyd  at  the  University  of  Michigan.  The  collide  library  takes  input 
from  DRACO  of  the  particles  grouped  by  cells  and  then  groups  the  particles  into  collision 
pairs,  calculates  the  probability  of  collision  for  each  pair,  performs  the  collisions  and  then 
returns  the  particles  to  DRACO  with  post-collision  velocities.  Upgrading  DRACO  for 
using  the  collide  library  requires  that  DRACO  has  the  ability  to  group  particles  by  cells. 
DRACO  does  not  typically  group  particles  by  cells  but  instead  simply  tracks  the  particles 
position  and  velocity  and  deposits  the  particles  charge  onto  the  nodes  of  the  Cartesian 
mesh.  An  algorithm  was  written  to  group  particles  by  cells  and  using  the  collide  library’s 
data  structure  to  perform  DSMC  collisions. 

A  test  simulation  was  run  using  the  collide  library  after  being  integrated  into 
DRACO  and  the  code  was  shown  to  produce  the  proper  number  of  collisions.  Further 
analysis  of  the  new  DSMC  code  is  preformed  as  part  of  this  thesis  and  the  results  are 
compared  to  MCC  and  the  differences  are  discussed.  The  results  of  this  analysis  are 
discussed  in  Chapter  5. 

3.3:  Recombination  Collisions 

The  COLISEUM  code  has  been  upgraded  to  be  able  to  perform  recombination 
collisions  in  plasma  plumes.  A  recombination  collision  occurs  when  an  ion  and  an 
electron  collide  to  form  a  neutral  particle  and  typically  only  occur  in  high  densities.  There 
are  two  types  of  recombination  collisions  of  interest  for  use  in  COLISEUM:  two-body 
and  three-body  collisions.  A  two-body  recombination  collision  is  a  direct  collision 
between  an  ion  and  an  electron  and  results  in  a  neutral  particle  and  a  photon.4 

Xe+  +e~  =>  Xe  +  hv  (3.7) 

The  procedure  for  performing  a  recombination  collision  first  checks  to  see  that  an 
elastic  collision  between  an  ion  and  an  electron  has  occurred.  For  each  collision  between 
an  ion  and  an  electron  the  probability  of  a  recombination  collision  is  calculated  and 
compared  to  a  random  number  between  zero  and  one.  If  the  random  number  is  less  than 
the  probability  then  a  recombination  collision  is  preformed.4 
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(3.10) 

(3.11) 


In  these  equations  the  input  variables  are  the  parameters  from  Bird’s  VHS 
collision  model  where  Tref  is  the  reference  gas  temperature,  dref  is  the  impact  parameter, 
mr  is  the  reduced  mass  of  the  collision,  £a  is  the  reaction  activation  energy  which  is  zero 

for  recombination,  £c  is  the  collision  energy  between  the  two  particles,  k  is  the 

Boltzmann  constant,  T  is  a  gamma  function  and  ar2  and  b,o  are  the  parameters  for  the 
modified  Arrhenius  rate  coefficient.  The  velocity  of  the  resulting  neutral  particle  is  a 
function  of  the  collision  particles  and  their  masses.4 
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(3.12) 


A  three-body  collision  occurs  when  an  ion  collides  with  an  electron  and  a 
secondary  electron  is  present  and  the  result  is  a  neutral  and  an  electron. 

Xe+  +e~  +e~  =>  Xe  +  e~  (3.13) 


The  probability  of  a  three-body  collision  is  calculated  in  a  similar  way  as  a  two- 
body  collision  with  a  few  subtle  differences.  A  three-body  collision  starts  with  an  elastic 
collision  between  an  ion  and  an  electron  and  a  randomly  paired  secondary  electron.  The 
probability  for  a  three-body  collision  uses  the  plasma  number  density,  ne.4 
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The  resulting  neutral  particle  has  a  velocity  equal  to  the  center-of-mass  velocity, 
ucm ,  and  secondary  electron  velocity,  uae . 
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(3.17) 
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The  input  parameters  used  to  test  these  recombination  probabilities  are  listed  in 
the  table  below  for  a  collision  between  a  Xenon  ion  and  an  electron. 


Table  3.1:  Input  Parameters  for  Recombination  Collisions 


CO 

0 

a 

2.3e-16 

b 

-0.5 

Tref 

300 

dref 

8.0e-10 

After  testing  the  added  recombination  code  to  COLISEUM  it  has  been  found  that 
the  probability  for  a  recombination  collision  is  very  small.  For  a  plasma  density  on  the 
order  of  1018  the  probability  of  a  recombination  collision  is  on  the  order  of  10"8  for  a  two- 
body  collision  and  10'10  for  a  three-body  collision.4  Obviously  a  very  high  plasma  density 
is  required  to  get  enough  recombination  collisions  to  have  any  effect  on  the  simulation. 
The  simulations  in  this  thesis  do  not  have  a  plasma  density  high  enough  to  have  any 
impact  on  the  simulations  and  in  fact  it  is  unlikely  that  any  recombination  collisions 
would  occur  during  the  entire  simulation.  As  a  result,  recombination  collisions  were  not 
considered  further  for  a  topic  in  this  thesis. 

3.4:  HPHall  Source 

A  new  source  has  been  developed  using  data  collected  from  the  Hall  Thruster 
device  code  HPHall,  originally  developed  by  M.  Fife,  in  order  to  further  enhance  the 
DRACO  code.  It  was  concluded  in  a  paper  presented  at  the  42nd  Joint  Propulsion 
Conference  that  the  source  has  a  very  significant  impact  on  the  results  of  PIC 
simulations.  Thus,  the  best  way  to  improve  the  results  of  DRACO  simulations  is  to  have 
the  best  source  representation  possible.  The  source  presented  in  that  paper, 
FLUX_R_VZ_VR,  is  based  on  experimental  laser  induced  fluoresce  (LIF)  data.  Although 
this  source  was  the  best  representation  of  a  Hall  Thruster  at  the  time,  the  results  did  not 
match  the  results  of  the  vacuum  chamber  experiment.  Complete  LIF  data  in  3-dimensions 
would  be  ideal  because  that  would  give  a  complete  description  of  the  velocity  distribution 
function  (VDF).  However,  the  200 W  source  is  based  on  data  taken  along  three  normal 
axis.  Thus,  it  is  not  possible  to  correlate  the  three  independent  components  and  the  source 
does  not  resolve  the  original  VDF.  The  data  from  HPHall  represents  correlated 
discretized  VDF  and  may  provide  a  more  accurate  representation  of  a  hall  thruster  exit 
plane. 
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3.4.1:  HPHall 


HPHall  is  a  hybrid-PIC  code  developed  to  study  the  physics  of  the  acceleration 
process  in  Hall  Thrusters.  HPHall  uses  quasi-one-dimensional  fluid  equations  for 
electrons  and  tracks  ions  and  neutral  particles  with  a  Boltzmann  solver.14  An  HPHall 
simulation  was  run,  using  the  input  parameters  for  the  200  W  Busek  Hall  Effect  Thruster, 
for  use  in  DRACO.  The  results  for  the  HPHall  simulation  are  show  in  the  figures  below. 
The  plots  are  contour  plots  of  potential,  temperature,  ion  particle  energy  and  ion  number 
density.  The  HPHall  code  can  also  output  other  parameters  not  shown  in  the  figures 
below.20 


Figure  3.4:  Contour  Plots  of  HPHall  Results:  Potential  (Top  Left),  Temperature  (Top  Right),  Particle 
Energy  (Bottom  Left)  and  Particle  Number  Density  (Bottom  Right) 
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3.4.2:  HPHall  Source  Algorithm 


The  particles  crossing  the  exit  plane  of  the  thruster  at  steady  state  are  output  into  a 
file  at  the  end  of  a  HPHall  simulation  and  the  data  is  then  input  into  DRACO.  The 
particle  information  that  is  output  from  HPHall  is  the  z-location  away  from  the  exit  plane, 
the  radial  location,  the  z-velocity  and  the  radial  velocity.  A  small  sample  of  the  output  is 
shown  below: 

VARIABLES  =  z  r  vz  vr 
ZONE  T=HPHALL_xel 

0.0205325  0.0113228  3365.72  -1073.27 

0.0205275  0.0103901  2761.25  -847.351 

The  z-location  is  ignored  in  this  algorithm  and  it  is  assumed  that  the  particles  start 
at  the  exit  plane  of  the  thruster.  The  HPHall  data  is  one-dimensional  in  position  and  two 
dimensional  in  velocity  and  must  be  modified  for  use  in  DRACO  to  have  a  3-dimensional 
position  and  velocity.  The  particle  information  is  transferred  into  DRACO  using  a  series 
of  rotations  and  coordinate  frames.  The  source  centroid  and  normal  vector  are  obtained  in 
DRACO  and  are  used  to  establish  the  source  coordinate  frame  (S-Frame).  The  three  unit 
vectors  used  to  establish  the  S-frame  are  the  unit  normal  vector,  n ,  a  vector  from  the 
centroid  to  a  node  on  a  triangular  surface  element,  f , ,  and  the  vector  cross  product 
between  the  normal  vector  and  the  tangent  vector,  f2 .  The  unit  vectors  that  make  up  this 
coordinate  frame  written  in  the  inertial  frame  are: 


?1 =  TJ  +  TlyJ  +  r J 

+  (3.19) 

n  —  ni  +  «  ,  /  +  n,k 

x  y  J  z 


The  next  coordinate  frame  used  is  the  particle  coordinate  frame  (P-Frame)  at 
which  the  first  unit  vector  of  the  frame  points  to  the  position  of  the  particle  that  is 
introduced  with  the  radius  and  velocity  from  the  HPHall  data  and  a  random  angle  from 
the  source  frame.  This  random  angle  is  the  device  used  to  transform  the  one-dimensional 
position  and  two-dimensional  velocity  to  3-dimensional  position  and  velocity.  This 
rotation  is  a  simple  rotation  about  the  3-axis  and  can  be  preformed  using  the  following 
simple  rotation  matrix  for  an  Euler  angle  rotation  about  the  3-axis.17 
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Figure  3.5:  Example  Source  Geometry  with  S-Frame  and  P-Frame  unit  vectors 


The  rotation  matrix  from  the  S-Frame  to  the  inertial  frame  can  be  written  simply 
as  the  three  mutually  orthogonal  unit  vectors  of  the  S-Frame  written  in  the  inertial  frame. 


[IS]  = 


nx 


'lx 


n„ 


2z 

n. 


(3.21) 


The  rotation  matrix  from  the  P-Frame  to  the  inertial  frame  can  be  written  as  the 
product  of  the  rotation  matrix  from  the  P-Frame  to  the  S-Frame  and  the  rotation  matrix 
from  the  S-Frame  to  the  inertial  frame. 


fax  C0S  6  -  Z2x  sin  0)  fax  sin  6  +  Z2x  C0S  0) 

nx 

[IP]= 

fay  cos  0  -  Tly  sin &)  fay  sin  6  +  r2y  cos  o) 

ny 

faz  cos  6  -  t2z  sin  6)  faz  sin  d  +  t2z  cos  0) 

nz _ 

The  vector  to  the  particle  expressed  in  the  P-Frame  is  simply  the  radius,  r, 
obtained  from  the  HPHall  data  in  the  /),  direction.  The  p-vector  can  be  expressed  in  the 
inertial  frame  as  the  product  of  the  rotation  matrix  with  the  p-vector  written  in  the  P- 
Frame.  Additionally,  since  the  source  is  not  necessarily  at  the  origin,  the  centroid 
position,  x0 ,  must  be  added  to  the  position  vector  as  well. 

pp  =rpl+0p2+0p3  (3.23) 

p‘=[lP]pp+x0  (3.24) 

The  velocity  vector  expressed  in  the  P-Frame  is  written  as  vr  taken  from  the 
HPHall  data  in  the  px  direction  and  vz  taken  from  the  HPHall  data  in  the  p  ,  direction. 
Additionally,  after  observing  the  initial  results  of  the  source,  a  theta  component  was 
added  to  the  velocity  of  all  the  particles,  ve,  in  the  p2.  All  of  the  particles  receive  the 
same  theta  component  equal  to  the  thermal  velocity.  Again,  the  velocity  vector  written  in 
the  inertial  frame  is  the  product  of  the  velocity  vector  and  the  rotation  matrix. 


19 

Distribution  A:  Approved  for  public  release;  distribution  unlimited 


(3.25) 


Vp  =vrPi+vep2+vzpi 

v'  =[ip]vp  (3.26) 

The  new  HPHall  source  is  specified  in  the  coliseum  input  file  and  requires  the 
particle  specie  name,  mass  flow  rate,  theta  velocity  and  the  file  name  of  the  HPHall 
output  file.  The  source  can  be  specified  in  the  input  file  in  the  following  way: 

source_specify  [source  component]  HPHALL  [particle  specie]  [  m  ]  [ v0]  [file  name] 

ex.  source_specify  sourcel  HPHALL  XE+  le-5  224.1  part_ionl.dat 


20 

Distribution  A:  Approved  for  public  release;  distribution  unlimited 


Chapter  4:  Design  of  Simulation 


An  experiment  was  conducted  at  the  Air  Force  Research  Laboratory  at  Edwards 
Air  Force  Base  in  California  and  the  data  gathered  was  used  to  validate  the  DRACO 
code.  A  200W  Busek  Hall  Effect  Thruster  was  used  during  the  experiment  and  a  variety 
of  probes  were  used  to  collect  data.  The  DRACO  simulation  is  designed  to  recreate  the 
AFRL  experiment  so  that  the  results  can  be  compared  directly  with  the  data  collected. 

4.1 :  Air  Force  Research  Laboratory  Experiment 

The  AFRL  experiment  discussed  in  this  thesis  was  conducted  by  Nakles  et.  al.  at 
the  Electric  Propulsion  Facility  at  Edwards  Air  Force  Base  in  Chamber  6. 19  The  vacuum 
chamber  is  a  cylindrical  tank  with  a  diameter  of  1.8  meters  and  a  length  of  2.5  meters. 
The  pressure  of  the  tank  during  the  experiment  is  approximately  7x1  O'4  Pa  and  the 
vacuum  is  maintained  by  two  cryogenic  vacuum  pumps.12 

4.1.1:Busek  Hall  Thruster 

The  thruster  used  during  the  experiment  is  a  Busek  BHT-200  Hall  Effect  Thruster 
with  an  input  power  of  200  W  and  a  thrust  of  12.8  mN.9  The  thruster  has  a  mass  flow  rate 
of  Xenon  propellant  of  0.94  mg/s  and  a  specific  impulse  of  1390  s.9  The  thruster  has  a 
discharge  voltage  of  250  V  and  an  efficiency  of  approximately  43%.9 


Figure  4.1:  Photograph  of  Busek  BHT-200  Hall  Thruster9 


A  HyperMesh  version  of  the  200  W  thruster  was  created  at  the  Air  Force 
Research  Laboratory  for  use  in  this  thesis.  The  thruster  is  divided  into  four  components: 
thruster,  thruster  exit  plane,  cathode  and  cathode  exit  plane.  Ions  are  injected  into  the 
simulation  from  the  thruster  exit  plane  and  neutrals  can  be  injected  into  the  simulation 
from  the  thruster  exit  plane. 
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Figure  4.2:  Plot  of  Surface  Mesh  of  Busek  BHT-200  Hall  Thruster 


4.1.2:  Experiment  Probes 


A  series  of  probes  were  used  to  gather  the  data  in  the  AFRL  experiments 
including  a  Faraday  probe,  Langmuir  probe,  double  probe  and  ExB  probe.  All  of  the 
probes  were  mounted  on  the  same  translation  system  which  moves  in  the  radial  and 
angular  directions.  A  photograph  of  the  translation  setup  is  included  below  inside  AFRL 
Chamber  6  with  the  Double  probe  attached.  The  probe  sweeps  from  -90  degrees  to  90 
degrees  in  radial  increments  from  2-12  cm,  depending  on  the  probe,  to  60  cm.  Multiple 
sweeps  for  each  probe  were  performed  to  insure  accurate  data.19 


Translation 

System 


Cryogenic 
Vacuum  Pump 


Double 

Probe 


Busek  200W 
Hall  Thruster 


Figure  4.3:  Photograph  of  Probe  Translation  System  in  AFRL  Chamber  619 
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A  guarded  Faraday  probe  was  used  to  collect  current  density  data.  A  Faraday 
probe  is  biased  below  the  plasma  potential  to  repel  electrons  away  from  the  probe  and 
measure  only  ion  current  into  the  probe.  The  probe  used  in  the  experiments  was  biased  to 
-30  V  with  respect  to  the  chamber  ground.  The  inner  collector  of  the  probe  has  a  diameter 
of  0.312”  and  the  outer  shield  has  a  diameter  of  0.885”. 19 


A  Langmuir  probe  was  used  to  measure  the  plasma  potential  and  electron 
temperature  of  the  Hall  thruster  plasma  plume.  A  Langmuir  probe  works  by  inserting  an 
electrode  into  the  plasma  with  an  electric  potential.  The  probe  had  a  wire  that  was 
thoriated  tungsten  0.015”  in  diameter  and  a  length  of  0.5”.  Originally,  a  emissive  probe 
was  suppose  to  be  used  for  gathering  plasma  potential  and  electron  temperature  data  but 
the  heated  wire  produced  noisy  data.  Additionally,  a  double  Langmuir  probe  was  used  to 
collect  ion  density  data  during  the  experiment.  The  probe  is  comprised  of  two  tungsten 
rods,  2”  long  and  0.15”  in  diameter  spaced  0.5”  apart.19 
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The  final  probed  used  in  the  AFRL  experiments  of  the  200W  Busek  Hall  Thruster 
was  an  ExB  probe.  This  probe  uses  uniform  perpendicular  electric  and  magnetic  fields 
that  are  both  normal  to  the  particle  velocity  vector.  These  fields  are  adjusted  such  that 
there  is  no  net  force  on  a  particle  with  a  given  velocity.  As  a  result,  only  particles  of  the 
given  velocity  can  travel  into  the  probe  and  thus  the  probe  can  show  the  current  collected 
by  the  probe  for  each  velocity.  The  probe  used  in  the  AFRL  experiments  is  detailed  in 
Reference  7. 

4.2:  Simulation  Set-up 

Most  of  the  DRACO  simulations  run  in  this  thesis  have  the  same  mesh  as  well  as 
the  same  thruster  geometry.  The  thruster  exit  plane  is  centered  at  the  origin  of  the  domain 
and  the  x-axis  and  y-axis  span  from  -0.30  m  to  0.30  m  and  the  z-axis  spans  from 
-0.08  m  to  0.60  m.  The  cells  of  the  mesh  are  cubes  and  have  a  volume  of  one  cubic 
centimeter.  The  mesh  contains  60  cells  in  the  x  and  y  directions  and  68  cells  in  the  z 
direction.  A  plot  of  the  simulation  mesh  as  well  as  the  empty  domain  with  the  thruster 
geometry  is  shown  below.  The  only  simulations  with  a  different  mesh  are  the  test  cases 
involving  mesh  effects  and  these  meshes  have  a  stretched  mesh  with  a  varying  cell  size. 
The  reference  values  used  by  the  field  solver  are  derived  from  the  experimental  data.  The 
time  step  used  by  the  simulation  is  set  to  an  automatic  adjustment  algorithm  such  that  the 
fastest  moving  particles,  in  this  case  the  double  charged  ions,  do  not  move  further  than  a 
fraction  of  a  cell.  For  these  cases  the  automatic  time  step  is  set  such  that  particles  to  not 
move  further  than  25%  of  a  cell  during  a  single  iteration.  The  criterion  of  25%  of  a  cell 
was  determined  through  a  brief  sensitivity  analysis  where  the  automatic  time  step 
parameter  was  varied  until  the  simulation  results  were  independent  of  the  time  step. 


Figure  4.6:  Plot  of  Simulation  Mesh  (Left)  and  Domain  with  Thruster  Mesh  (Right) 
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The  simulations  are  set  to  run  to  a  steady  state  solution  based  on  the  number  of 
particles  entering  and  exiting  the  domain.  When  the  number  of  particles  entering  the 
domain  is  approximately  equal  to  the  number  exiting  the  domain  then  DRACO  assumes 
that  the  simulation  has  reached  steady  state.  After  steady  state  the  simulations  run  for  an 
additional  thousand  iterations  to  ensure  that  the  simulation  results  are  indeed  steady  state 
and  then  yet  another  additional  thousand  iterations  of  probe  sampling  and  averaging  of 
the  potential  field  data.  Sampling  for  a  thousand  iterations  after  steady  state  provides  a 
good  time  average  of  the  data  with  a  reasonable  noise  level. 

4.2.1 :  Simulation  Probes 

Two  probes  have  been  added  to  DRACO  to  validate  the  code  against 
experimental  probe  data.  A  current  density  probe  was  written  to  compare  Faraday  probe 
data  and  an  ExB  probe  was  written  to  compare  ExB  probe  data. 

The  DRACO  current  density  probe  first  checks  to  see  if  a  particle  is  past  the 
sampling  radius  and  then  calculates  the  angle  with  respect  to  the  thruster  exit.  The  charge 
of  the  particle  is  then  deposited  on  the  probe  in  the  corresponding  probe  angle.  The 
sampling  time  is  tracked  as  well  as  the  probe  area  and  used  to  calculate  the  current 
density  at  the  end  of  the  simulation.  The  sampling  angles  and  current  density  is  output 
into  a  file  at  the  end  of  the  simulation  and  then  plotted  and  compared  with  experimental 
data.  The  angles  are  output  in  degrees  and  the  current  density  is  output  in  amps  per 
square  meter. 

A  simple  ExB  probe  has  been  written  in  DRACO  to  sample  the  current  for  a 
range  of  particle  velocities.  The  probe’s  function  is  simple;  a  particle’s  position  is 
compared  with  the  position  and  angle  of  the  probe  to  see  if  the  particle  is  in  the  sample 
range.  If  the  particle  is  in  range  of  the  probe  then  the  particles  velocity  is  compared  to  the 
probe’s  view  angle.  If  the  probe  is  within  the  view  angle  then  the  particles  charge  is  noted 
in  the  probe  for  the  corresponding  velocity  of  the  particle.  The  sampling  time  is  tracked 
and  at  the  end  of  the  simulation  the  charge  is  divided  by  the  sampling  time  to  get  the 
current.  The  data  is  output  into  a  file  of  velocities  and  current  and  is  then  plotted  and 
compared  to  the  experimental  data.  Both  sets  of  data  are  later  post-processed  to  have 
velocities  in  units  of  electron  volts  and  currents  that  are  normalized  on  a  scale  from  zero 
to  one. 

4.3:  Test  Cases 

A  sensitivity  study  was  preformed  to  validate  the  DRACO  code  and  to  determine 
what  input  parameters  have  the  greatest  effect  on  the  results  of  the  simulation.  A  series  of 
test  cases  were  conducted  to  investigate  all  of  the  important  input  parameters.  The  first 
step  in  running  the  simulations  is  to  determine  the  reference  parameters  to  put  into  the 
input  files.  To  do  this  the  experimental  data  was  investigated  to  determine  the  best  first 
approximation  for  the  reference  parameters.  After  the  first  approximation  was  made  the 
data  was  put  into  Tecplot  and  the  reference  parameters  were  tweaked  until  the  simulation 
and  experimental  plots  fit  best.  Then  these  reference  parameters  were  used  for  the 
reminder  of  the  simulations. 
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4.3.1 :  Particle  Source 


It  was  discovered  in  a  previous  paper  presented  in  the  42nd  Joint  Propulsion 
Conference  that  the  source  model  has  a  major  impact  on  the  results  of  a  simulation.25  To 
improve  the  results  of  the  simulation  a  new  source  was  created  based  on  the  Hall  Thruster 
device  code  HPHall.  The  particle  source  is  the  first  input  parameter  investigated  because 
it  has  such  a  large  impact  on  the  simulation.  The  HPHall  source  was  compared  to  the 
FLUX_R_VZ_VR  source  which  was  the  previous  state-of-the-art  source  used  in  previous 
paper  presented  at  JPC.  These  sources  are  also  compared  to  the  MAXSTREAM  source 
which  is  a  very  basic  source  based  on  a  Maxwellian  velocity  distribution.  Three 
variations  of  the  HPHall  source  are  also  compared  to  determine  how  to  best  model  a  Hall 
Thruster  plume  using  the  new  source.  The  first  variation  tested  is  the  standard  HPHall 
source  that  samples  at  the  exit  plane  of  the  thruster.  The  next  variation  samples  at  the  end 
of  the  HPHall  domain  rather  than  the  exit  plane  of  the  thruster.  The  final  variation 
samples  at  the  exit  plane  of  the  thruster  but  also  loads  the  z-direction  electric  field  from 
HPHall  into  DRACO.  These  three  variations  are  compared  to  each  other  as  well  as 
experimental  data  to  determine  how  each  variation  affects  the  steady  state  plasma  plume 
results. 

4.3.2:  Field  Solver 

There  are  several  different  field  solving  methods  that  can  be  implemented  in 
DRACO.  No  non-linear  solvers  will  be  used  in  this  thesis  but  the  parameters  of  the 
Boltzmann  solver  will  be  examined  to  determine  the  importance  of  solving  the  field  in 
these  simulations.  Simulations  will  be  run  without  solving  the  electric  field  as  well  as 
simulations  using  the  Boltzmann  solver  with  the  constant  and  polytropic  temperature 
models.  The  results  of  the  field  solving  simulations  will  be  compared  to  experimental 
results  to  determine  how  the  results  different  and  which  field  solving  method  produces 
the  best  results. 

4.3.3:  Collision  Model 

There  are  two  primary  methods  for  performing  particle  collisions  in  DRACO: 
Monte  Carlo  Collisions  and  Direct  Simulation  Monte  Carlo.  Typically,  the  MCC  methods 
are  faster  than  DSMC  methods  but  do  not  conserve  momentum.  There  are  also  two 
collision  types  investigated  in  this  study:  Charge  Exchange  and  Variable  Hard  Sphere.  A 
series  of  simulations  are  conducted  to  investigate  the  effects  of  particle  collisions  on  the 
results.  First  simulations  are  conducted  without  any  particle  collisions  at  all.  Second 
particles  are  conducted  with  MCC  and  CEX  collisions  using  a  projected  neutral 
background  gas.  The  third  test  case  uses  MCC  with  CEX  and  VHS  collisions  while 
actually  tracking  the  neutral  particles.  Finally,  DSMC  is  used  with  CEX  and  VHS 
collisions. 
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4.3.4 :  Mesh-Effects 


The  recently  added  mesh-stretching  DRACO  code  is  tested  in  this  sensitivity 
study  to  determine  when  stretching  can  be  applied  and  how  the  results  are  affected  by  the 
mesh  stretching.  A  linear  stretching  scheme  is  applied  in  all  three  directions  such  that  the 
mesh  is  finest  at  the  exit  plane  of  the  thruster  and  coarsest  at  the  boundary  of  the  domain. 
The  results  are  compared  to  the  uniform  mesh  as  well  as  experimental  data. 


Table  4.1:  List  of  Test  Cases 


Test  Case 

Description 

la 

HPHall  Source,  MCC  CEX  Collisions 

lb 

FLUX_R_VZ_VR  Source,  MCC  CEX  Collisions 

lc 

MAXSTREAM  Source,  MCC  CEX  Collisions 

Id 

HPHall  Modified  Source,  MCC  CEX  Collisions 

le 

HPHall  Z-Direction  Electric  Field  Loading  Source,  MCC  CEX  Collisions 

2a 

HPHall  Source,  No  Field  Solver,  MCC  CEX  Collisions 

2b 

HPHall  Source,  Boltzmann  Const,  MCC  CEX  Collisions 

2c 

HPHall  Source,  Boltzmann  Poly,  MCC  CEX  Collisions 

3a 

HPHall  Source,  No  Collisions 

3b 

HPHall  Source,  MCC  CEX  Collisions 

3c 

HPHall  Source,  MCC  CEX  &  VHS  Collisions 

3d 

HPHall  Source,  DSMC  CEX  &  VHS  Collisions 

4a 

HPHall  Source,  MCC  CEX  Collisions  with  Linear  Stretched  Mesh 

4b 

HPHall  Electric  Field  Loading  Source,  MCC  CEX  Collisions,  Linear  Stretched  Mesh 
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Chapter  5:  Simulation  Results 

The  results  of  the  many  simulations  conducted  as  part  of  this  thesis  are  displayed 
and  discussed  in  this  chapter.  This  includes  contour  plots  and  line  plots  of  the  plasma 
density,  polytropic  temperature  and  plasma  potential  as  well  as  probe  output  for  ExB  and 
Faraday  probes.  All  of  the  simulations  use  the  same  reference  parameters:  NeRef  = 
7.67xl017  #/m3,  teref  =  5  eV,  phiref  =  26  V,  and  gamma  =  1.2. 

5.1 :  DRACO  Source  Study 

The  first,  and  perhaps  the  most  important,  study  performed  in  this  thesis  examines 
the  affect  of  the  particle  source  used  in  DRACO  on  the  results  of  the  simulation.  There 
are  three  sources  being  considered  in  this  thesis:  MAXSTREAM,  FLUX_R_VZ_VR,  and 
the  newly  developed  HPHall  source.  Additionally,  multiple  variations  of  the  HPHall 
source  must  be  examined  to  determine  how  to  best  use  the  source. 


5.1.1:  MAXSTREAM,  FLUX  &  HPHall 


z,  m  z,  m 

(c)  FLUX_R_VZ_VR  (d)  HPHall 

Figure  5.1:  Contour  Plots  of  Plasma  Potential  for  Experiment  and  Sources 


28 

Distribution  A:  Approved  for  public  release;  distribution  unlimited 


The  plot  of  the  experimental  data  shows  a  higher  potential  than  any  of  the 
DRACO  simulations  directly  at  the  exit  plane  of  the  thruster.  For  radii  between  10  cm 
and  30  cm  the  experimental  data  appears  to  agree  within  2-3  V  for  all  three  sources.  In 
order  to  make  a  more  precise  analysis  of  the  data  the  contour  plots  are  converted  into  line 
plots.  The  angles  from  the  center  of  the  exit  plane  of  the  thruster  are  calculated  and  the 
potential  is  plotted  as  a  function  of  the  angle. 


Figure  5.2:  Plots  of  Plasma  Potential  as  a  Function  of  Angle  for  Experiment  and  Sources 


At  a  radius  of  6  cm  away  from  the  exit  plane  of  the  thruster  the  peak  value  of  the 
experimental  data  matches  the  HPHall  and  FLUX  sources  however  the  peak  is  much 
wider  than  the  simulations.  The  MAXSTREAM  source  peaks  at  a  much  lower  value  than 
the  experimental  data.  For  angles  higher  than  35  degrees  the  experimental  data  is  closest 
to  the  MAXSTREAM  source  and  in  between  the  FLUX  and  HPHall  source.  At  a  radius 
of  12  cm  away  from  the  exit  plane  the  experimental  data  near  the  center  of  the  thruster 
more  closely  matches  the  values  of  the  simulations.  The  peak  is  much  narrower  and  the 
curve  for  the  FLUX  and  HPHall  source  almost  exactly  matches  from  5  to  20  degrees.  For 
angles  higher  than  20  degrees  the  curve  matches  best  for  the  MAXSTREAM  source  and 
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values  in  between  the  FLUX  and  HPHall  source.  At  further  radii,  24  cm  and  30  cm,  the 
peak  values  for  the  experimental  data  closely  matches  the  MAXSTREAM  and  HPHall 
sources  but  the  FLUX  source  peaks  several  volts  higher  than  the  experimental  data. 
However,  for  angles  greater  than  30  degrees,  the  FLUX  source  produces  results  that  most 
closely  match  the  experimental  results. 


z,  m 

(a)  Experiment 


0.3 

z,  m 

(b)  MAXSTREAM 


(c)  FLUX_R_VZ_VR 

Figure  5.3:  Contour  Plots  of  Polytropic  Temperature  for  Experiment  and  Sources 


It  can  be  easily  seen  in  the  above  contour  plots  of  the  polytropic  temperature  that 
the  agreement  is  not  as  strong  as  that  of  the  plasma  potential.  The  shape  of  the  plume 
appears  similar  for  the  MAXSTREAM  and  HPHall  cases  but  the  values  are  all  higher 
than  that  of  the  experiment.  Where  as  the  FLUX  source  also  has  higher  values  than  that 
of  the  experiment  and  the  shape  of  the  plume  differs  as  well.  These  differences  can  be 
seen  more  clearly  in  line  plots  as  a  function  of  the  angle  from  the  center  of  the  thruster 
(Figure  5.4).  Very  near  the  thruster  exit  plane  the  polytropic  temperature  agrees  pretty 
well  for  the  experiment  and  the  simulations.  The  results  further  away  from  the  thruster 
exit  plane  produce  a  similar  curve  for  the  MAXSTREAM  and  HPHall  source  but  are  off 
by  a  scalar  factor. 
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Figure  5.4:  Plots  of  Polytropic  Temperature  as  a  Function  of  Angle  for  Experiment  and  Sources 


At  first  glance  it  appears  that  the  cause  of  the  disagreement  in  the  polytropic 
temperature  is  due  to  the  reference  values  used  for  solving  the  potential.  In  an  attempt  to 
correct  this  disagreement,  the  reference  temperature  was  varied  but  the  result  was  that  the 
polytropic  temperature  results  improved  and  the  plasma  potential  results  degraded.  This 
is  likely  due  to  the  fact  that  a  polytropic  temperature  model  assumes  a  constant  specific 
heat  ratio  (Equation  5.1).  If  the  specific  heat  ratio  is  calculated  using  the  experimental 
data  and  then  plotted  on  a  contour  plot  (Figure  5.5)  it  can  be  clearly  seen  that  the  specific 
heat  ratio  is  not  constant.  Near  the  exit  plane  of  the  thruster  the  value  of  the  specific  heat 
ratio  is  approximately  1.6  where  as  further  away  the  value  is  closer  to  1.2.  The  reference 
specific  heat  ratio  was  chosen  to  correspond  to  the  value  seen  in  the  low  density  charge 
exchange  wings. 


T  = 


f  V~l 

n 


\noj 


(5.1) 
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Figure  5.5:  Plot  of  Specific  Heat  Ratio,  y,  Calculated  from  Experimental  Data 
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(c)  FLUX_R_VZ_VR  (d)  HPHall 

Figure  5.6:  Contour  Plots  of  Ion  Number  Density  for  Experiment  and  Sources 
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Contour  plots  of  the  ion  number  density  (Figure  5.6)  show  a  good  agreement 
between  the  experimental  data  and  the  simulations  in  main  plume  for  angles  less  than  30 
degrees,  especially  for  HPHall  simulations.  It  should  be  noted  that  the  probe  used  for  the 
AFRL  experiments  for  the  particle  number  density  can’t  collect  data  closer  than  12  cm 
away  from  the  thruster  exit.  Again,  the  results  for  the  FLUX  simulations  show  a  narrower 
beam  than  that  of  the  experiment  and  other  simulations. 


(a)  r  =  12  cm 


Figure  5.7:  Plots  of  Ion  Number  Density  as  a  Function  of  Angle  for  Experiment  and  Sources 


The  line  plots  of  the  particle  number  density  as  a  function  of  angle  (Figure  5.7) 
show  a  similar  agreement.  At  a  radius  of  12  cm  the  curve  for  the  HPHall  simulations 
almost  lies  exactly  on  top  of  experimental  results  up  to  an  angle  of  70  degrees.  The 
results  for  all  three  simulations  peak  at  approximately  the  same  values  as  the  experiment 
for  all  radii  but  the  main  peak  of  the  FLUX  simulations  is  very  narrow  and  for  angles 
greater  than  30  degrees  the  density  is  basically  uniform.  The  results  for  higher  angles,  24 
cm  and  30  cm,  still  show  a  similar  peak  value  for  MAXSTREAM  and  HPHall  but  the 
charge  exchange  wings  level  off  at  a  value  an  order  of  magnitude  higher  then  that  of  the 
experimental  data.  The  results  of  plasma  potential,  polytropic  temperature  and  number 
density  show  that  these  plasma  properties  can  be  recreated  in  DRACO  in  the  main  beam 
using  any  of  the  three  sources. 
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Figure  5.8:  Plots  of  Faraday  Probe  Results  Compared  to  Experimental  Data 


The  Faraday  probe  analysis  clearly  shows  that  the  HPHall  source  closely  agrees 
with  the  experimental  data  over  the  entire  range  of  angles  and  radii.  The  entire  two  curves 
almost  overlap  except  in  the  direct  center  of  the  plume  where  the  HPHall  source  peaks 
slightly  higher  than  the  experimental  results.  The  next  closest  fit  is  the  MAXSTREAM 
source  which  has  a  very  similar  trend  as  the  experimental  data  but  the  charge  exchange 
wings  level  off  at  a  lower  current  density.  The  main  beam  of  the  MAXSTREAM  source 
also  nearly  overlaps  from  angles  of  0  to  approximately  40  degrees.  The  poorest 
agreement  to  the  experimental  results  is  the  FLUX  source  which  has  a  much  narrower 
main  beam  and  peaks  at  a  value  higher  than  the  experimental  data.  The  change  exchange 
wings  are  much  wider  and  level  off  at  a  current  density  lower  than  that  of  the 
experimental  data.  These  results  would  tend  to  indicate  that  the  HPHall  source  is  the  best 
of  the  three  sources  but  further  analysis  is  still  required. 

Next  the  ExB  probe  data  will  be  analyzed  to  determine  if  the  particles  are  being 
produce  with  the  correct  velocities  in  the  correct  quantities.  The  experimental  data 
provided  is  extensive  and  each  data  set  analyzed  provides  slightly  different  results.  The 
experimental  results  presented  in  this  paper  are  believed  to  be  results  that  have  been 
published. 
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(a)  MAXSTREAM 


(b)  FLUX_R_VZ_VR 


(c)  HPHall 


Figure  5.9:  Plots  of  ExB  Probe  Results  at  an  Angle  of  0  Degrees  Away  from  Thruster 


The  ExB  probe  results  above  show  that  the  highest  current  agrees  relatively  well 
with  the  experimental  results  for  the  FLUX  and  HPHall  sources.  Both  peaks  are  very 
slender  and  occur  between  150  eV  and  250  eV.  Although  this  appears  to  be  a  wide  range 
the  experimental  data  provided  shows  a  similar  range  in  sampling.  It  can  be  clearly  seen 
that  the  MAXSTREAM  source  produces  particles  that  move  too  fast  as  well  as  a  wider 
distribution  in  the  velocity  of  the  particles.  Additionally,  the  FLUX  and  MAXSTREAM 
sources  only  have  one  central  peak  and  no  sub-peaks  where  as  the  HPHall  source  has  a 
central  peak  and  a  sub-peak  that  represents  the  population  of  Xe++  ions.  This  sub-peak 
also  shows  up  in  the  experimental  data  and  can  be  clearly  seen  in  these  plots.  In  other 
data  sets  provided  this  sub-peak  is  better  defined  and  in  others  it  is  less  defined  but  this 
peak  appears  at  the  expected  location,  twice  the  energy  of  the  main  peak.  The 
experimental  results  also  show  a  third  peak  at  triple  the  energy  of  the  main  peak  and  this 
represents  the  population  of  Xe+++  ions,  which  are  not  considered  in  the  simulations.  It 
is  clear  from  these  results  that  the  HPHall  source  produces  results  that  best  match  the 
experimental  results  and  the  next  best  results  come  from  the  FLUX  source.  The  FLUX 
source  results  can  also  be  easily  modified  by  changing  the  values  in  the  input  file  to  more 
closely  match  the  results  of  this  particular  set  of  experimental  data.  However,  only  one 
peak  would  still  exist.  The  ExB  results  for  MAXSTREAM  have  poor  agreement. 
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(a)  MAXSTREAM 


(b)  FLUX_R_VZ_VR  (c)  HPHall 

Figure  5.10:  Plots  of  ExB  Probe  Results  at  an  Angle  of  10  Degrees  Away  from  Thruster 


(a)  MAXSTREAM 


(b)  FLUX_R_VZ_VR  (c)  HPHall 

Figure  5.11:  Plots  of  ExB  Probe  Results  at  an  Angle  of  20  Degrees  Away  from  Thruster 
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(a)  MAXSTREAM 


(b)  FLUX_R_VZ_VR  (c)  HPHall 

Figure  5.12:  Plots  of  ExB  Probe  Results  at  an  Angle  of  30  Degrees  Away  from  Thruster 


The  ExB  probe  analysis  for  increasing  sampling  angles  (Figure  5.10  -  Figure 
5.12)  shows  that  the  MAXSTREAM  and  FLUX  simulations  remain  approximately  the 
same  and  the  HPHall  source  results  diverge  from  the  experimental  results  for  increasing 
sample  angle.  The  HPHall  source  at  a  sampling  angle  of  10  degrees  away  from  the 
thruster  exit  has  a  similar  agreement  as  the  results  for  0  degrees.  Where  as  at  20  and  30 
degrees  the  main  peak  becomes  wider  and  the  secondary  peak  begins  to  become  less 
defined  as  it  melds  into  the  primary  peak.  Also,  at  30  degrees  the  peak  value  occurs  at  a 
lower  energy  than  in  previous  simulations.  The  main  particle  population  peaks  at  a  low 
energy  range,  less  than  100  eV. 

This  analysis  shows  that  the  newly  developed  HPHall  source  produces  results  that 
most  closely  model  a  200  W  Busek  Hall  Effect  Thruster  in  a  vacuum  chamber.  The  best 
agreement  can  be  seen  from  the  Faraday  probe  results  where  the  experimental  and 
simulation  curves  for  the  HPHall  source  almost  overlap  for  the  entire  curve.  The  ExB 
probe  results  show  that  the  HPHall  source  can  capture  additional  information,  such  as 
sub-peaks,  that  the  other  sources  can  not.  However,  further  analysis  is  still  required  to 
determine  how  to  best  utilize  the  new  source. 
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5.1.2:  HPHall  Source  Variations 


A  further  investigation  of  the  HPHall  source  has  been  conducted  to  determine 
how  to  best  use  the  source.  Multiple  sets  of  experimental  data  were  collected  at  AFRL 
and  three  of  the  data  sets  are  shown  in  Figure  5.13.  It  can  be  seen  in  the  data  set  displayed 
in  the  previous  section,  Figure  5.13a,  that  some  of  the  data  agrees  very  well  with  the 
DRACO  simulations  where  as  the  other  two  data  sets  do  not  agree  as  well.  In  Figure 
5.13b  the  main  peak  is  located  at  a  slightly  higher  energy  and  the  two  sub-peaks  occur  at 
a  much  higher  energy  and  are  more  defined  than  the  other  data  sets.  Figure  5.13c  shows 
that  the  main  peak  occurs  at  almost  double  the  energy  of  the  main  peak  of  the  simulation 
and  the  sub-peaks  are  not  nearly  as  defined. 


Energy,  eV 

(a) 


Energy,  eV  Energy,  eV 


(b)  (c) 

Figure  5.13:  Plots  of  ExB  Probe  Results  for  HPHall  Source  and  Various  Experimental  Data  Sets 


In  an  attempt  to  improve  the  results  produced  by  the  HPHall  source  two  variations 
of  the  source  have  been  investigated.  Upon  further  investigation  of  the  HPHall  simulation 
output  a  large  z-direction  electric  field  (Figure  5.14),  on  the  order  of  105  V/m,  can  be 
seen  just  outside  the  exit  plane  of  the  thruster  (z  =  20  mm).  This  is  hypothesized  to  be  the 
reason  why  the  particles  are  not  being  accelerated  as  much  as  the  experimental  data.  The 
two  possible  solutions  to  this  problem  are:  load  the  particles  into  the  simulations  at  the 
end  of  the  HPHall  simulation  or  load  the  electric  field  from  HPHall  into  DRACO. 
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Figure  5.14:  Contour  Plots  of  HPHall  Z-Direction  Electric  Field  Output  (Left)  and  Overlaid 
Cartesian  Mesh  past  Thruster  Exit  Plane  (Right) 


The  first  method  tested  was  sampling  the  particles  at  the  outer  boundary  of  the 
HPHall  simulations,  as  shown  in  the  blue  particles  in  Figure  5.15.  A  new  source  mesh  is 
required  for  use  with  this  method  which  would  match  the  boundary  for  which  the 
particles  are  being  sampled  (Figure  5.16).  This  is  a  major  disadvantage  because  the  actual 
geometry  of  the  thruster  can  not  be  modeled  in  the  experiment  and  it  is  unknown  what  the 
potential  of  the  surface  should  be  set  to.  Additionally,  particles  will  be  injected  into  the 
simulation  at  a  point  in  the  middle  of  the  domain  instead  of  directly  outside  of  the  thruster 
exit  plane.  This  also  means  that  any  physics  that  occurs  between  the  exit  plane  of  the 
thruster  and  the  sampling  geometry  will  be  lost  if  it  is  not  handled  in  HPHall,  including 
charge  exchange  collisions. 
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Figure  5.15:  Plot  of  Particle  Sampling  at  the  end  of  the  HPHall  Domain 
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Figure  5.16:  Plot  of  Mesh  used  with  Modified  HPHall  Sampling 


The  second  method  tested  loads  the  z-direction  electric  field  produced  by  the 
HPHall  code  into  DRACO.  This  is  done  by  use  of  a  MATLAB  script  created  to 
interpolate  the  electric  field  onto  the  DRACO  Cartesian  mesh.  The  interpolated  mesh  can 
be  seen  in  the  figure  below.  It  can  obviously  be  seen  in  Figure  5.17  that  the  resolution  of 
the  current  DRACO  mesh  does  not  provide  a  very  good  resolution  of  the  electric  field. 
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(b)  Y-Z  Slice  Zoom 
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Figure  5.17:  Contour  Plots  of  Z-Direction  Electric  Field  Input  into  DRACO 
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Although  the  resolution  of  the  electric  field  is  not  well  captured  in  the  uniform 
mesh  it  is  still  used  in  this  study  for  testing.  The  plots  of  the  potential  for  the  experiment, 
standard  HPHall  source,  HPHall  source  that  samples  at  the  outer  boundary  and  the 
HPHall  source  that  loads  the  z-direction  electric  field  are  shown  in  the  figure  below. 
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(c)  HPHall  Outer  Boundary  Sampling  (d)  HPHall  Loaded  Electric  Field 

Figure  5.18:  Contour  Plots  of  Plasma  Potential  for  Experiment  and  HPHall  Source  Methods 

It  can  be  seen  in  the  contour  plots  in  Figure  5.18  that  the  plasma  potential  is 
affected  by  the  variation  of  the  HPHall  source  used.  The  plot  of  the  HPHall  source  that 
samples  at  the  outer  boundary  of  the  HPHall  simulation  clearly  shows  a  gap  between 
where  the  particles  are  injected  into  the  simulation  and  the  thruster  exit  plane.  The 
general  shape  of  the  plume  is  also  different  from  the  normal  HPHall  source.  The  HPHall 
source  that  loads  the  electric  field  from  HPHall  produces  a  main  beam  in  the  near  plume 
region  that  appears  to  more  closely  resemble  the  experimental  data.  However,  due  to  the 
poor  resolution  of  the  loaded  electric  field  there  exists  an  asymmetry  in  the  plume  which 
does  not  show  up  in  the  other  cases  or  the  experimental  data.  This  could  be  fixed  by  using 
a  finer  mesh  resolution  or  a  non-uniform  mesh  and  will  be  further  investigated  in  a  later 
section  of  this  thesis.  To  more  closely  investigate  the  differences  between  the  sources  the 
contour  plots  are  again  converted  into  line  lots  as  a  function  of  angle  from  the  thruster. 
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Figure  5.19:  Line  Plots  of  Potential  as  a  function  of  Angle  from  Thruster  for  HPHall  Source  Methods 


The  line  plots  of  plasma  potential  in  Figure  5.19  show  that  at  a  radius  of  6  cm  all 
of  the  sources  have  approximately  the  same  agreement  with  the  experimental  data. 
However,  at  radii  from  12  cm  to  30  cm  the  HPHall  source  that  utilizes  a  loaded  z- 
direction  electric  field  produces  results  that  most  closely  resemble  the  experimental  data. 
At  12  cm,  the  curve  for  the  experimental  data  and  the  field  loading  HPHall  source 
overlap  for  the  majority  of  angles,  except  for  angles  less  than  10  degrees  or  greater  than 
70  degrees.  The  agreement  for  the  central  beam  begins  to  diverge  at  24  cm  and  30  cm  and 
the  curves  no  longer  overlap  until  angles  of  approximately  25  degrees.  The  agreement  for 
the  standard  HPHall  source  and  the  modified  sampling  HPHall  source  have 
approximately  the  same  agreement  and  do  not  match  as  well  as  the  field  loading  HPHall 
source. 


The  contour  plots  of  ion  density  shown  in  Figure  5.20  have  similar  differences 
between  the  three  HPHall  source  variations  as  the  plasma  potential  plots.  Again  we  see  in 
the  other  boundary  sampling  simulation  that  the  particles  are  entering  the  domain  away 
from  the  thruster  exit  plane  and  the  electric  field  loading  source  has  asymmetries  in  the 
plume. 
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(c)  HPHall  Outer  Boundary  Sampling 


(d)  HPHall  Loaded  Electric  Field 


Figure  5.20:  Contour  Plots  of  Ion  Number  Density  for  Experiment  and  HPHall  Source  Methods 
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(b)  r  =  12  cm 


(c)  r  =  24  cm 


(c)  r  =  24  cm 


Figure  5.21:  Line  Plots  of  Ion  Number  Density  as  a  function  of  Angle  from  Thruster  for  HPHall 

Source  Methods 
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The  line  plots  of  ion  number  density  as  a  function  of  angle  (Figure  5.21)  show 
that  at  a  radius  of  12  cm  away  from  the  thruster  the  results  are  best  for  the  standard 
HPHall  source  and  the  electric  field  loading  HPHall  source.  However,  at  24  cm  and  30 
cm  the  results  for  angles  from  10  degrees  to  40  degrees  are  best  for  the  field  loading 
HPHall  source.  At  angles  above  50  degrees  for  a  radius  of  24  cm  and  30  cm  the  ion 
number  density  produced  by  all  three  variations  of  the  HPHall  source  do  not  closely 
match  the  experimental  data.  The  plasma  potential  and  ion  number  density  plots  tend  to 
indicated  that  the  electric  field  loading  HPHall  source  produces  the  best  results.  Next, 
probe  results  will  be  displayed  and  discussed. 


Figure  5.22:  Plots  of  Faraday  Probe  Results  for  HPHall  Source  Methods 


The  Faraday  probe  results  shown  in  Figure  5.22  show  that  the  best  fit  to  the 
experimental  data  is  still  the  standard  HPHall  source.  The  modified  sampling  HPHall 
source  also  agrees  well  with  the  Faraday  probe  data  but  the  electric  field  loading  HPHall 
sources  does  not  agree  as  well  with  the  experimental  data.  The  electric  field  loading 
source  produces  current  densities  lower  than  the  experimental  data  for  angles  greater  than 
20  degrees  but  still  has  a  curve  with  a  similar  shape.  This  could  be  due  to  the 
asymmetries  from  the  poor  resolution  of  the  electric  field  loaded  into  DRACO. 
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(b)  HPHall  Outer  Boundary  Sampling  (c)  HPHall  Loaded  Electric  Field 

Figure  5.23:  Plots  of  ExB  Probe  Results  at  an  Angle  of  10  Degrees  Away  from  Thruster  for  HPHall 

Source  Methods 


It  can  be  clearly  seen  in  Figure  5.23  that  the  ExB  probe  results  for  the  outer 
boundary  sampling  HPHall  source  and  the  electric  field  loading  HPHall  source  peak  at  a 
higher  energy  than  the  standard  HPHall  source.  Additionally,  the  main  peaks  of  these  two 
methods  are  wider  than  the  main  peak  of  the  standard  HPHall  source  and  the 
experimental  results.  The  cause  of  the  wider  peak  in  the  electric  field  loading  HPHall 
source  could  be  due  to  the  poor  resolution  of  the  electric  field  loaded  into  DRACO  which 
causes  the  particles  to  not  all  be  accelerated  to  the  correct  velocity.  In  addition,  as 
discussed  above,  the  energy  of  the  main  peak  varies  depending  on  the  experimental  data 
set  used.  Thus,  all  three  sources  have  a  main  peak  in  the  general  energy  range  expected  to 
be  produced  by  this  thruster. 

This  analysis  and  the  analysis  in  the  previous  section  shows  that  for  a  uniform 
mesh  the  standard  HPHall  source  produces  results  that  most  closely  resemble  the 
experiment  conducted  at  AFRL.  As  a  result  the  simulations  in  the  future  sections  will  use 
the  standard  HPHall  source.  The  mesh  study  will  include  an  analysis  using  the  electric 
field  loading  HPHall  source  as  well  because  a  finer  mesh  near  the  exit  plane  of  the 
thruster  will  likely  improve  the  results.  The  results  produced  by  the  HPHall  source  most 
closely  resemble  the  experimental  data  taken  from  a  200  W  Busek  Hall  Effect  Thruster 
and  will  greatly  enhance  the  capability  of  the  DRACO  code. 
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5.2:  Field  Solver  Study 

The  next  study  conducted  as  part  of  this  thesis  is  a  sensitivity  analysis  that 
investigates  the  affects  of  the  field  solver  used  during  the  simulation.  Three  variations  of 
the  field  solver  are  simulated:  Boltzmann  solver  with  a  polytropic  temperature  model, 
Boltzmann  solver  with  a  constant  temperature  model  and  no  field  solver  at  all. 
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(b)  Boltzmann  Polytropic  Temperature  Model  (c)  Boltzmann  Constant  Temperature  Model 

Figure  5.24:  Contour  Plots  of  Plasma  Potential  for  Experiment  and  Field  Solving  Methods 


The  contour  plots  of  the  plasma  potential  for  the  experimental  results  and  the 
simulations  with  the  various  field  solving  methods  are  shown  above.  Obviously  the 
simulation  with  no  field  solver  does  not  require  a  contour  plot  because  all  of  the  potential 
values  are  zero.  The  center  of  the  plume  near  the  thruster,  less  than  20  cm,  has  similar 
results  for  both  the  poly  tropic  and  constant  temperature  model  simulations.  However,  the 
constant  temperature  model  shows  a  much  smaller  and  narrower  plume  with  positive 
plasma  potential  values  and  most  of  the  values  are  zero  or  less.  At  first  glance  it  appears 
that  the  polytropic  temperature  model  produces  results  that  have  significantly  better 
agreement  to  the  experimental  data  than  that  of  the  constant  temperature  model.  To  verify 
this,  line  plots  are  again  created  to  more  directly  compare  the  results  (Figure  5.25). 
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Figure  5.25:  Line  Plots  of  Plasma  Potential  for  Experiment  and  Field  Solving  Methods 


It  can  be  seen  in  the  above  plots  that  the  polytropic  and  constant  temperature 
models  produce  very  similar  results  in  the  central  beam  very  near  the  thruster  exit.  At  a 
radius  of  6  cm  the  curves  for  the  polytropic  and  constant  temperature  models  overlap 
until  approximately  20  degrees  where  the  two  curves  diverge.  At  a  radius  of  12  cm  the 
curves  for  poly  tropic  and  constant  temperature  models  both  peak  at  approximately  15  V 
with  the  experimental  data.  However,  now  the  two  curves  diverge  for  angles  greater  than 
5  degrees.  The  poly  tropic  temperature  model  curve  closely  resembles  the  experimental 
data  for  a  radius  greater  than  12  cm  and  differs  only  by  at  most  3  Y  throughout  the  range 
of  angles.  For  a  radius  of  24  cm  or  greater  the  constant  temperature  model  does  not  match 
the  polytropic  temperature  model  or  experimental  data  at  all  and  the  values  for  angles 
greater  than  approximately  15  degrees  are  negative.  All  of  the  constant  temperature 
model  curves  tend  to  be  much  more  linear  than  the  polytropic  temperature  model  and  the 
experimental  data.  As  noted  in  the  contour  plots  the  curve  for  no  field  solver  is  zero  over 
the  entire  range  of  angles  and  radii. 
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Figure  5.26:  Contour  Plots  of  Polytropic  Temperature  for  Experiment  and  Field  Solving  Methods 


Figure  5.27:  Line  Plots  of  Poly  tropic  Temperature  for  Experiment  and  Field  Solving  Methods 
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The  contour  plots  for  the  polytropic  and  constant  temperature  models  are  almost 
identical  where  as  the  results  without  a  field  solver  differ.  The  center  of  the  plume  has  a 
longer  range  of  higher  temperature  and  away  from  the  center  of  the  plume  the 
temperatures  are  lower  than  the  other  simulations.  The  line  plots  show  that  at  6  cm  the 
polytropic  and  constant  temperature  models  produce  almost  identical  temperature  curves 
where  as  the  no  field  simulation  produces  a  curve  of  similar  shape  but  approximately 
0.5eV  higher  than  the  other  simulation  curves.  At  12  cm,  the  no  field  case  produces  a 
curve  that  no  longer  has  the  same  shape  as  the  other  simulation  curves.  The  peak  of  the 
curve  differs  from  the  other  simulation  curves  more  than  the  previous  radius  and  the  drop 
in  the  curve  is  more  dramatic  than  the  other  simulations.  Again,  at  24  cm  and  30  cm  the 
peak  is  more  pronounced  and  then  for  angles  greater  than  30  degrees  the  temperature  is 
lower  than  the  other  simulations.  For  angles  greater  than  75  degrees  the  no  field  case 
produces  results  that  are  also  lower  than  the  experimental  data.  In  general,  the  results  for 
the  polytropic  and  constant  temperature  models  produce  similar  curves  as  the 
experimental  results  but  at  a  higher  temperature,  where  as  the  no  field  results  do  not 
closely  match  the  experimental  results  at  all. 
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(b)  Boltzmann  Polytropic  Temperature  Model 


(c)  Boltzmann  Constant  Temperature  Model 


(d)  No  Field  Solver 


Figure  5.28:  Contour  Plots  of  Ion  Number  Density  for  Experiment  and  Field  Solving  Methods 
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The  contour  plots  of  ion  number  density  above  again  show  that  the  polytropic  and 
constant  temperature  models  show  similar  results.  The  no  field  simulation  results  show  a 
narrow  plume  with  a  higher  density  in  the  central  plume  than  the  other  simulations  and 
the  experimental  results.  The  line  plots  below  show  that  the  results  for  all  three 
simulations  agree  relatively  well  at  a  radius  of  12  cm  for  angles  less  than  70  degrees.  The 
simulation  with  no  field  solver  produces  results  with  higher  densities  than  the  other 
simulations  and  the  experiment  and  the  polytropic  and  constant  temperature  models 
produce  nearly  identical  results.  For  a  radius  of  24  cm  and  30  cm  the  polytropic  and 
constant  temperature  models  produce  similar  results  as  the  experiment  for  angles  less 
than  30  degrees  but  higher  densities  in  the  charge  exchange  wings  by  approximately  an 
order  of  magnitude.  The  simulations  with  no  field  solver  produce  results  in  the  center  of 
the  beam  which  are  higher  than  the  other  simulations  and  experiment  and  than  does  not 
produce  a  curve  that  matches  any  of  the  other  results.  The  results  of  the  polytropic 
temperature  and  ion  number  density  show  that  the  polytropic  and  constant  temperature 
models  both  produce  similar  results  and  that  the  results  with  no  field  solver  do  not  match 
well  to  the  other  simulations  or  the  experiment.  Next  the  probe  results  will  be  analyzed 
and  discussed. 


Figure  5.29:  Line  Plots  of  Ion  Number  Density  for  Experiment  and  Field  Solving  Methods 
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Figure  5.30:  Faraday  Probe  Results  for  Experiment  and  Field  Solving  Methods 


The  Faraday  probe  results  (Figure  5.30)  show  that  the  polytropic  and  constant 
temperature  models  are  nearly  identical  over  the  entire  curve  and  both  closely  match  the 
experimental  results.  At  a  radius  of  10  cm  the  curve  with  no  field  solving  peaks  at 
approximately  the  same  value  as  the  other  simulations  but  does  not  match  any  of  the 
other  curves.  For  angles  greater  than  60  degrees  the  agreement  between  the  simulation 
with  no  field  solver  and  the  other  curves  differ  by  over  an  order  of  magnitude  or  more. 

The  ExB  probe  results  (Figure  5.31  -  Figure  5.34)  again  show  that  the  polytropic 
and  constant  temperature  models  produce  similar  results.  The  constant  temperature  model 
peaks  at  a  value  closer  to  the  experimental  results  than  the  polytropic  model.  The  case 
with  no  field  solver  shows  the  main  peak  with  a  wider  distribution  than  the  other  cases  or 
the  experimental  data  and  peaks  at  a  lower  energy.  The  secondary  peak  is  also  less 
defined  for  the  case  with  no  field  solver.  At  higher  angles,  20  -  30  degrees,  the 
simulation  results  for  all  three  cases  become  more  similar,  although  the  case  with  no  field 
solver  has  more  noise  than  the  other  cases  despite  the  fact  that  the  sampling  time  was  the 
same  for  all  of  the  simulations. 
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(a)  Boltzmann  Polytropic  Temperature  Model 


(b)  Boltzmann  Constant  Temperature  Model  (c)  No  Field  Solver 

Figure  5.31:  ExB  Probe  Results  at  an  Angle  of  0  degrees  for  Experiment  and  Field  Solving  Methods 


(a)  Boltzmann  Polytropic  Temperature  Model 


(b)  Boltzmann  Constant  Temperature  Model 


(c)  No  Field  Solver 


Figure  5.32:  ExB  Probe  Results  at  an  Angle  of  10  degrees  for  Experiment  and  Field  Solving  Methods 
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(a)  Boltzmann  Polytropic  Temperature  Model 


(b)  Boltzmann  Constant  Temperature  Model  (c)  No  Field  Solver 

Figure  5.33:  ExB  Probe  Results  at  an  Angle  of  20  degrees  for  Experiment  and  Field  Solving  Methods 


(a)  Boltzmann  Polytropic  Temperature  Model 


(b)  Boltzmann  Constant  Temperature  Model  (c)  No  Field  Solver 

Figure  5.34:  ExB  Probe  Results  at  an  Angle  of  30  degrees  for  Experiment  and  Field  Solving  Methods 
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The  results  of  the  field  study  show  that  the  polytropic  and  constant  temperature 
models  produce  similar  results  for  the  majority  of  the  plasma  properties  and  probe  results. 
The  exception  is  the  plasma  potential  which  differs  greatly  between  the  two  cases  and  the 
results  for  the  polytropic  temperature  model  more  closely  match  the  experimental  results 
than  the  constant  temperature  model.  Additionally,  the  above  results  show  that  the  field 
solver  is  a  necessary  part  of  the  simulation  to  produce  satisfactory  results.  The  case  with 
no  field  solver  does  not  produce  results  that  agree  well  with  the  experimental  results.  This 
case  also  produces  noisy  results  which  may  be  due  to  the  random  number  generator  used 
for  the  simulations  in  this  thesis.  The  lack  of  field  solver  means  that  the  particles  move  in 
a  straight  line  after  entering  the  simulation  and  thus  the  random  number  generator  has 
more  influence  over  this  case.  As  a  result  of  this  analysis,  the  Boltzmann  solver  with  the 
polytropic  temperature  model  is  used  for  the  remainder  of  the  simulations  in  this  thesis. 

5.3:  Particle  Collision  Study 

The  next  study  conducted  as  part  of  this  thesis  is  a  sensitivity  analysis  of  particle 
collisions  in  DRACO  simulations.  As  discussed  in  Reference  18,  particle  collisions  can 
have  a  major  impact  on  the  results  of  a  simulation  using  a  Hall  Thruster  in  a  vacuum 
chamber.  This  study  will  examine  four  cases: 

1 .  No  Particle  Collisions 

2.  MCC  model  with  a  projected  background  neutral  density  and  CEX  collisions  only 

3.  MCC  model  with  neutral  particle  tracking  and  both  CEX  and  VHS  collisions 

4.  DSMC  model  with  neutral  particle  tracking  and  both  CEX  and  VHS  collisions 

Neutral  particles  can  be  handled  two  ways  in  DRACO:  projected  background  density 
and  tracking  the  actual  neutral  particles.  For  the  projected  neutral  background  density  an 
analytical  model  for  the  neutrals  is  generated  for  the  source  at  the  start  of  the  simulation 
and  a  background  density  is  specified  for  the  neutrals.  As  a  result  you  get  an 
approximation  of  the  steady  state  neutral  density  for  use  during  the  entire  simulation  and 
you  do  not  need  to  track  neutral  particles  during  the  simulation.  This  method,  however, 
only  works  for  CEX  collisions  because  VHS  collisions  require  momentum  exchange 
between  two  particles. 

Simulations  that  track  neutral  particles  can  use  both  MCC  and  DSMC  collision 
methods  and  can  track  both  CEX  and  VHS  collisions.  However,  simulations  that  track 
neutral  particles  are  more  computationally  demanding  and  the  neutral  particles  loaded  at 
the  start  of  the  simulation  are  removed  as  the  simulation  runs.  Neutral  particles  that  reach 
the  boundary  of  the  domain  are  removed  from  the  simulation  as  well  as  neutrals  that  are 
part  of  a  CEX  collision.  As  a  result  the  neutral  density  of  the  simulation  is  not  constant  as 
it  is  with  the  neutral  projection  method. 

The  collision  simulations  will  all  use  the  Boltzmann  field  solver  with  the  polytropic 
temperature  model  and  the  standard  HPHall  source.  The  background  neutral  density  used 
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in  the  simulations  is  2xl017  particles/m3  (7xl0~4  Pa)  and  neutrals  that  come  out  from  the 
source  due  to  inefficiencies  are  accounted  for  as  10%  of  the  source  ion  mass  flow  rate. 


z,  m 


(a)  Experiment 


z,  m  z,  m 


(d)  MCC  CEX-VHS  (e)  DSMC  CEX-VHS 

Figure  5.35:  Contour  Plots  of  Plasma  Potential  for  Experiment  and  Collision  Study  Simulations 
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The  contour  plots  of  plasma  potential  for  the  collision  study  (Figure  5.35)  show 
that  all  of  the  simulations  are  nearly  identical  in  the  main  beam  near  the  thruster  exit.  It  is 
also  noteworthy  that  the  contour  plots  for  both  of  the  MCC  cases  have  similar  plumes  as 
well  as  the  cases  for  no  collisions  and  the  DSMC  case.  The  no  collision  case  and  the 
DSMC  case  both  have  a  noisy  plume  with  similar  values  throughout  the  plume.  The  line 
plots  of  plasma  potential  (Figure  5.36)  show  that  for  a  radius  of  6  cm  and  12  cm,  all  of 
the  simulations  produce  the  same  curve  for  angles  less  than  20  degrees.  At  a  radius  of  24 
cm  and  30  cm  the  curves  match  for  angles  less  than  10  degrees.  Although  the  curves 
don’t  match  at  higher  angles  the  differences  between  the  curves  is  only  1  -  3  V.  Again, 
the  curves  for  no  collisions  and  the  DSMC  test  case  match  closely  throughout  the  entire 
range  of  angles,  especially  at  lower  radii.  The  agreement  between  the  simulation  results 
and  the  experimental  data  is  approximately  the  same  for  all  of  the  cases. 


Figure  5.36:  Line  Plots  of  Plasma  Potential  for  Experiment  and  Collision  Study  Simulations 
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(a)  Experiment 
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(b)  No  Collisions  (c)  MCC  CEX  Only 


(d)  MCC  CEX-YHS  (e)  DSMC  CEX-YHS 


Figure  5.37:  Contour  Plots  of  Polytropic  Temperature  for  Experiment  and  Collision  Study 

Simulations 
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Figure  5.38:  Line  Plots  of  Polytropic  Temperature  for  Experiment  and  Collision  Study  Simulations 


The  results  for  polytropic  temperature  (Figure  5.37  -  Figure  5.38)  show  a  similar 
agreement  and  trends  as  the  results  for  the  plasma  potential.  Again  it  would  seam  that  the 
results  near  the  exit  of  the  thruster  in  the  main  beam  are  not  greatly  affected  by  the 
difference  in  collision  models  but  the  differences  are  increased  as  the  distance  and  angle 
from  the  thruster  is  increased.  It  is  again  noteworthy  to  point  out  the  similarities  between 
the  no  collision  and  the  DSMC  simulation. 
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Figure  5.39:  Contour  Plots  of  Ion  Number  Density  for  Experiment  and  Collision  Study  Simulations 
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Figure  5.40:  Line  Plots  of  Ion  Number  Density  for  Experiment  and  Collision  Study  Simulations 


The  results  for  ion  number  density  (Figure  5.39  -  Figure  5.40)  again  show  similar 
agreement  and  trends  as  the  plasma  potential  and  polytropic  temperature.  Again  the 
results  for  no  collisions  and  the  DSMC  simulation  are  most  similar.  The  most  plausible 
explanation  as  to  why  the  no  collision  and  DSMC  simulations  produce  similar  results  is 
that  at  steady  state  the  DSMC  simulation  is  not  producing  enough  collisions.  This  is 
likely  because  the  neutrals  are  being  removed  from  the  simulation  and  the  steady  state 
neutral  density  is  too  low.  Contour  plots  of  the  neutral  density  at  the  end  of  the 
simulations  (Figure  5.41)  show  that  the  neutral  density  in  the  main  beam  is  more  than  an 
order  of  magnitude  lower  than  at  the  start  of  the  simulation  for  the  DSMC  case.  This  is 
also  true  for  the  MCC  simulation  that  track  both  CEX  and  VHS  collisions  but  not  to  the 
extent  as  the  DSMC  simulation.  If  the  number  of  collisions  per  time  step  is  plotted  as  a 
function  of  time  step  (Figure  5.42)  it  can  be  seen  that  the  DSMC  simulation  produces  the 
most  collisions  of  all  the  collision  simulations  but  the  number  of  collisions  decreases  as 
the  simulation  runs  rather  than  leveling  off. 
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(b)  MCC  CEX-VHS  (c)  DSMC  CEX-VHS 

Figure  5.41:  Contour  Plots  of  Neutral  Density  for  Experiment  and  Collision  Study  Simulations 
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Figure  5.42:  Plot  of  Number  of  Collision  per  Time  Step  as  a  function  of  Time  Step 
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Figure  5.43:  Plots  of  Faraday  Probe  Results  for  Experiment  and  Mesh  Simulations 


The  Faraday  probe  results  (Figure  5.43)  for  all  simulation  cases  are  basically 
identical  for  angles  less  than  50  degrees  over  the  entire  range  of  radii.  The  main 
differences  between  the  cases  occur  in  the  charge  exchange  wings,  approximately  50  -  90 
degrees.  In  this  region  the  only  simulation  that  produces  a  curve  that  closely  matches  the 
experimental  data  is  the  MCC  using  the  projected  background  neutral  density.  The  other 
curves  basically  continue  the  same  trend  as  the  main  beam  and  decay  exponentially  over 
the  entire  range  of  angles. 

The  ExB  probe  results  (Figure  5.44  -  Figure  5.47)  for  the  collision  study 
simulations  show  no  major  differences  in  the  results  between  the  cases  over  the  entire 
range  of  angles.  At  higher  angles  there  are  some  small  differences  in  the  main  peak  but 
the  peak  still  occurs  at  the  same  value. 
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(a)  No  Collisions  (b)  MCC  CEX  Only 


(c)  MCC  CEX-YHS  (d)  DSMC  CEX-YHS 


Figure  5.44:  Plots  of  ExB  Probe  Results  for  Experiment  and  Collision  Simulations  at  0  degrees 


(a)  No  Collisions 


(c)  MCC  CEX-YHS  (d)  DSMC  CEX-VHS 

Figure  5.45:  Plots  of  ExB  Probe  Results  for  Experiment  and  Collision  Simulations  at  10  degrees 
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(a)  No  Collisions 


(b)  MCC  CEX  Only 


(c)  MCC  CEX-YHS 


(d)  DSMC  CEX-YHS 


Figure  5.46:  Plots  of  ExB  Probe  Results  for  Experiment  and  Collision  Simulations  at  20  degrees 
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(c)  MCC  CEX-YHS  (d)  DSMC  CEX-VHS 

Figure  5.47:  Plots  of  ExB  Probe  Results  for  Experiment  and  Collision  Simulations  at  30  degrees 
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According  to  the  plasma  potential,  polytropic  temperature  and  ion  number  density 
results  the  collision  simulations  are  approximately  the  same,  especially  in  the  main  beam. 
However,  the  plots  of  neutral  density  (Figure  5.41)  and  collisions  as  a  function  of  time 
step  (Figure  5.42)  show  that  the  method  for  tracking  neutrals  can  have  major  implications 
on  the  steady  state  simulation  results.  When  particle  tracking  is  performed,  instead  of 
using  a  projected  neutral  density,  the  number  of  neutral  particles  in  the  simulation  is 
reduced  as  the  simulation  runs.  The  Faraday  probe  results  show  the  most  noticeable 
differences  between  the  collision  study  cases.  In  the  charge  exchange  wings  the 
simulations  that  track  neutral  particles  do  not  have  a  defined  wing  where  as  the 
simulation  using  a  projected  neutral  density  with  change  exchange  collision  tracking  only 
produces  a  curve  very  similar  to  that  of  the  experimental  data. 

The  results  from  this  study  show  that  for  a  vacuum  chamber  simulation  with  a 
high  neutral  density,  tracking  particle  collisions  is  important  to  produce  accurate  results 
in  the  charge  exchange  wings  of  a  Hall  thruster  plume.  The  results  also  show  that  the  best 
results  can  be  achieved  using  a  projected  background  neutral  density  instead  of  tracking 
neutral  particles.  Not  only  is  this  the  most  accurate  of  the  collision  models  but  it  is  also 
the  least  computationally  demanding  because  the  number  of  particles  in  the  simulation  is 
greatly  reduced.  As  a  result,  MCC  with  a  neutral  projection  is  the  collision  model  that  is 
used  during  the  remainder  of  this  thesis  and  this  was  also  the  collision  model  used  during 
the  source  and  field  study  cases. 
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5.4:  Mesh  Study 


The  final  study  conducted  as  part  of  this  thesis  is  a  mesh  sensitivity  analysis.  A 
uniform  and  non-uniform  mesh  are  tested  and  compared  to  determine  the  effects  of  the 
mesh  on  the  results  of  the  simulation.  The  uniform  mesh  has  60  cells  in  the  x  and  y 
directions  and  68  cells  in  the  z  direction  with  cubic  centimeter  cells.  The  non-uniform 
mesh  uses  a  linear  cell  size  stretching  routine  that  generates  a  mesh  with  the  finest  cell 
sizes  near  the  exit  plane  of  the  thruster.  The  x  and  y  direction  cells  are  2.5  cm  near  the 
outside  of  the  domain  and  0.5  cm  in  the  center  of  the  domain.  The  z  direction  mesh  is 
uniform  with  cell  length  of  0.5  cm  until  the  exit  plane  of  the  thruster  where  a  linear 
stretched  mesh  starts  with  cell  length  0.5  cm  and  at  the  end  of  the  domain  the  cell  length 
is  3.5  cm.  The  non-uniform  mesh  has  40  cells  in  the  x  and  y  directions  and  46  cells  in  the 
z  direction. 


(a)  Uniform  Mesh  (b)  Non-Uniform  Stretched  Mesh 

Figure  5.48:  Plots  of  Uniform  and  Non-Uniform  Meshes 


The  non-uniform  mesh  will  be  tested  for  the  standard  HPHall  source  and  the  z- 
direction  electric  field  loading  HPHall  source.  The  purpose  of  testing  both  sources  is  to 
both  demonstrate  the  capability  of  using  a  non-uniform  mesh  in  a  simulation  and  also 
show  how  a  non-uniform  mesh  can  be  used  to  refine  the  electric  field  loaded  into 
DRACO  from  HPHall. 

5.4.1 :  HPHall  Standard  Source 

The  first  non-uniform  mesh  simulation  that  was  run  used  the  standard  HPHall 
source.  The  potential  contour  plots  (Figure  5.49)  for  the  uniform  and  non-uniform  mesh 
show  that  the  results  for  both  simulations  are  almost  identical.  The  line  plots  of  potential 
(Figure  5.50)  also  show  that  the  curves  for  the  uniform  and  stretched  mesh  basically 
overlap  over  the  entire  range  of  angles.  The  biggest  differences  between  the  two  curves 
occur  for  angles  less  than  10  degrees  where  the  stretched  mesh  produces  a  curve  that 
more  closely  resembles  the  experimental  data.  This  is  likely  because  the  stretched  mesh  is 
finest  in  lower  angles  near  the  exit  of  the  thruster. 
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Figure  5.49:  Contour  Plots  of  Plasma  Potential  for  Uniform  and  Non-Uniform  Meshes 


Figure  5.50:  Line  Plots  of  Plasma  Potential  for  Experiment  and  Mesh  Simulations 
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(a)  Uniform  Mesh  (b)  Non-Uniform  Mesh 

Figure  5.51:  Contour  Plots  of  Polytropic  Temperature  for  Uniform  and  Non-Uniform  Meshes 


(a)  r  =  06  cm  (b)  r  =  12  cm 


(c)  r  =  24  cm  (d)  r  =  30  cm 

Figure  5.52:  Line  Plots  of  Poly  tropic  Temperature  for  Experiment  and  Mesh  Simulations 
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Figure  5.53:  Contour  Plots  of  Ion  Number  Density  for  Uniform  and  Non-Uniform  Meshes 
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Figure  5.54:  Line  Plots  of  Ion  Number  Density  for  Experiment  and  Mesh  Simulations 
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The  results  for  the  polytropic  temperature  and  ion  number  density  (Figure  5.51  - 
Figure  5.54)  show  that  the  results  for  the  stretched  mesh  are  almost  identical  to  the  results 
of  the  uniform  mesh.  Again,  the  biggest  difference  between  the  two  curves  occurs  at 
angles  below  10  degrees  and  all  the  differences  are  generally  small.  The  Faraday  probe 
results  (Figure  5.55)  show  the  same  agreement  as  the  previous  results  and  the  curve  for 
the  stretched  mesh  for  angles  less  than  10  degrees  is  a  better  fit  to  the  experimental 
results  than  the  uniform  mesh. 


Figure  5.55:  Plots  of  Faraday  Probe  Results  for  Experiment  and  Mesh  Simulations 


The  ExB  probe  results  (Figure  5.56)  show  that  the  uniform  and  non-uniform 
meshes  produce  similar  results  for  a  range  of  angles  from  0  to  20  degrees.  At  an  angle  of 
0  and  10  degrees  the  main  peak  is  identical  but  the  sub-peak  is  slightly  different  for  the 
non-uniform  mesh.  At  an  angle  of  20  degrees  the  main  peak  has  a  slightly  different  shape 
than  the  uniform  mesh.  The  non-uniform  mesh  results  show  a  sub-peak  at  an  energy 
approximately  50  eV  less  than  the  main  peak  where  as  the  uniform  mesh  simply  shows  a 
wide  main  peak.  Despite  these  differences,  the  ExB  probe  results  for  both  cases  are  very 
similar. 
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(a)  Uniform  Mesh,  0  Degrees  (b)  Non-Uniform  Mesh,  0  Degrees 


(c)  Uniform  Mesh,  10  Degrees  (d)  Non-Uniform  Mesh,  10  Degrees 


(e)  Uniform  Mesh,  20  Degrees  (f)  Non-Uniform  Mesh,  20  Degrees 

Figure  5.56:  Plots  of  ExB  Probe  Results  for  Experiment  and  Mesh  Simulations 


The  results  from  this  study  show  that  a  stretched  mesh  can  be  successfully 
implemented  and  produce  results  as  accurate  as  a  uniform  mesh.  This  functionally  can  be 
used  to  improve  results  with  a  finer  mesh  in  high  plasma  density  regions  as  well  as 
improve  computational  performance  due  to  the  need  for  fewer  cells.  The  uniform  mesh 
has  244,800  cells  where  as  the  non-uniform  mesh  has  73,600  cells,  thus  the  non-uniform 
mesh  has  30%  as  many  cells  as  the  uniform  mesh.  Obviously,  for  larger  domains  this 
functionality  can  be  extremely  useful. 
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5.4.2:  HPHall  Electric  Field  Loading  Source 


The  next  simulation  conducted  as  part  of  this  study  uses  the  z-direction  electric 
field  loading  HPHall  source.  As  discussed  in  the  source  study,  mesh  refinement  is 
required  when  loading  the  electric  field  from  HPHall  into  DRACO  to  improve  the  results. 
The  non-uniform  mesh  designed  for  use  in  this  study  has  the  finest  mesh  at  the  exit  plane 
of  the  thruster  where  the  electric  field  is  loaded.  Contour  plots  of  the  z-direction  electric 
field  for  the  uniform  and  non-uniform  meshes  are  shown  in  Figure  5.57.  It  is  obvious  that 
the  electric  field  loaded  in  the  non-uniform  mesh  is  considerably  more  detailed  than  the 
uniform  mesh.  It  can  also  be  seen  in  Figure  5.14  that  the  non-uniform  mesh  produces  an 
electric  field  that  is  more  similar  to  the  HPHall  electric  field. 
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Figure  5.57:  Contour  Plots  of  Z-Direction  Electric  Field  for  Uniform  and  Non-Uniform  Meshes 


The  contour  and  line  plots  of  potential,  polytropic  temperature  and  ion  number 
density  (Figure  5.58  -  Figure  5.63)  are  shown  below  and  like  the  results  using  the 
standard  HPHall  source  the  results  between  the  two  meshes  are  very  similar. 
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Figure  5.58:  Contour  Plots  of  Plasma  Potential  for  Uniform  and  Non-Uniform  Meshes 


Figure  5.59:  Line  Plots  of  Plasma  Potential  for  Experiment  and  Mesh  Simulations 
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(a)  Uniform  Mesh  (b)  Non-Uniform  Mesh 

Figure  5.60:  Contour  Plots  of  Polytropic  Temperature  for  Uniform  and  Non-Uniform  Meshes 


Figure  5.61:  Line  Plots  of  Polytropic  Temperature  for  Experiment  and  Mesh  Simulations 
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Figure  5.62:  Contour  Plots  of  Ion  Number  Density  for  Uniform  and  Non-Uniform  Meshes 
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Figure  5.63:  Line  Plots  of  Ion  Number  Density  for  Experiment  and  Mesh  Simulations 
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The  results  for  potential,  polytropic  temperature  and  ion  number  density  show  that 
the  uniform  and  non-uniform  meshes  produce  very  similar  results.  The  main  difference 
between  the  two  cases  is  that  the  uniform  mesh  has  more  noise  in  the  results  as  well  as  a 
more  defined  asymmetry  in  the  contour  plots.  This  is  likely  due  to  the  better  refined 
electric  field  loaded  into  DRACO  from  HPHall  in  the  non-uniform  mesh  case.  It  can  be 
assumed  that  further  mesh  refinement  would  likely  completely  remove  the  asymmetry 
from  the  simulation.  This  is  intuitively  obvious  because  the  electric  field  has  a  major 
impact  on  the  motion  of  the  particles  and  the  asymmetries  in  the  electric  field  would 
cause  asymmetries  in  the  simulation  results. 


Figure  5.64:  Plots  of  Faraday  Probe  Results  for  Experiment  and  Mesh  Simulations 


The  Faraday  probe  results  in  Figure  5.64  show  that  the  uniform  and  stretched 
mesh  produce  similar  results,  although  both  simulations  do  not  have  a  good  correlation  to 
the  experimental  results.  The  ExB  probe  results  in  Figure  5.65  show  a  noticeable 
improvement  between  the  uniform  and  non-uniform  meshes. 
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(a)  Uniform  Mesh,  0  Degrees  (b)  Non-Uniform  Mesh,  0  Degrees 


(c)  Uniform  Mesh,  10  Degrees  (d)  Non-Uniform  Mesh,  10  Degrees 


(e)  Uniform  Mesh,  20  Degrees  (f)  Non-Uniform  Mesh,  20  Degrees 

Figure  5.65:  Plots  of  ExB  Probe  Results  for  Experiment  and  Mesh  Simulations 


The  main  difference  between  the  ExB  probe  results  for  the  uniform  and  non- 
uniform  meshes  is  in  the  main  peak.  The  main  peak  for  the  non-uniform  mesh  is 
narrower,  especially  at  angles  of  0  and  10  degrees.  Additionally  the  sub-peak  is  better 
defined  for  the  stretched  mesh  case  as  well.  This  is  also  related  to  the  electric  field  mesh 
refinement  because  the  particles  are  receiving  a  more  equal  acceleration,  thus  the  main 
peak  is  narrower  because  the  particles  have  less  of  a  range  of  velocity. 
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The  results  in  this  test  case  show  that  if  a  stretched  mesh  is  correctly  implemented 
it  can  produce  results  that  are  as  accurate  as  a  uniform  mesh  and  even  more  so  in  areas  of 
most  refinement.  If  a  stretched  mesh  is  used  when  loading  an  electric  field  into  a 
simulation  it  can  improve  the  overall  results  of  the  simulation,  reducing  noise  and 
asymmetries  produced  by  a  poorly  refined  mesh  around  the  loaded  electric  field. 
Additionally,  the  reduction  of  cells  in  a  stretched  mesh  can  decrease  the  computational 
demands  of  a  simulation.  In  general  the  z-direction  electric  field  loading  HPHall  source 
does  not  produce  better  results  than  the  standard  HPHall  source.  As  stated  in  the  source 
study,  the  ExB  probe  experimental  data  varies  and  the  main  peak  for  the  electric  field 
loading  case  matches  some  data  sets  better  than  the  standard  HPHall  source.  For  the 
stretched  mesh  the  electric  field  loading  HPHall  source  also  produces  a  narrower  main 
peak  and  further  refinement  may  produce  a  peak  more  similar  to  the  experimental  data. 
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Chapter  6:  Conclusions 


This  thesis  discusses  the  recent  developments  to  the  electric  propulsion  plume 
code  DRACO  as  well  as  a  validation  and  sensitivity  analysis  of  the  code  using  data  from 
an  AFRL  experiment  using  a  Busek  200  W  Hall  Thruster.  DRACO  is  a  PIC  code  that 
models  particles  kinematically  while  using  finite  differences  schemes  to  solver  the 
electric  potential  and  field.  As  part  of  this  thesis,  a  new  particle  source  for  the  DRACO 
code,  based  on  the  Hall  Thruster  Device  code  HPHall,  was  developed.  Additionally,  the 
DRACO  code  was  developed  to  non-uniform  mesh  functionality,  the  capability  to 
perform  DSMC  collisions  and  recombination  collisions  together  with  Lubos  Brieda  and 
Alex  Barrie. 

The  DRACO  code  has  been  recently  modified  to  improve  simulation  results, 
functionality  and  performance.  A  particle  source  has  been  added  that  uses  the  Hall 
Thruster  device  code  HPHall  as  input  for  a  source  to  model  Hall  Thrusters.  This  source 
imports  2-dimensional  positions  and  velocities  of  neutral,  single  and  double  charged  ions 
from  HPHall  and  loads  them  into  DRACO  3-dimensionaltly  around  the  exit  plane  of  the 
thruster.  Multiple  variations  of  this  source  have  been  created  including  one  that  samples 
at  the  end  of  the  HPHall  simulation  and  one  that  loads  a  z-direction  electric  field  into 
DRACO  as  well  as  the  particles.  The  code  is  now  capable  of  using  a  non-uniform  mesh 
that  uses  any  combination  of  uniform,  linear  and  exponential  stretching  schemes  in  any  of 
the  three  directions.  A  stretched  mesh  can  be  used  to  refine  simulation  results  in  certain 
areas,  such  as  the  exit  of  a  thruster,  or  improve  performance  by  reducing  the  number  of 
cells  in  a  mesh.  Finally,  DRACO  now  has  the  capability  of  using  a  DSMC  collision 
scheme  as  well  as  performing  recombination  collisions. 

A  sensitivity  analysis  of  the  newly  upgraded  DRACO  code  was  performed  to  test 
the  new  functionalities  of  the  code  as  well  as  validate  the  code  using  experimental  data 
gathered  at  AFRL  using  a  Busek  200  W  Hall  Thruster.  A  simulation  was  created  that 
attempts  to  numerically  recreate  the  AFRL  experiment  and  the  validation  is  performed  by 
comparing  the  plasma  potential,  polytropic  temperature,  ion  number  density  of  the 
thruster  plume  as  well  as  Faraday  and  ExB  probe  results.  The  study  compares  the  newly 
developed  HPHall  source  with  older  source  models  and  also  compares  the  variations  of 
the  HPHall  source.  The  field  solver  and  collision  model  used  is  also  compared  to 
determine  how  to  achieve  the  best  results  using  the  DRACO  code.  Finally,  both  uniform 
and  non-uniform  meshes  are  tested  to  determine  if  a  non-uniform  mesh  can  be  properly 
implemented  to  improve  simulation  results  and  performance. 

The  validation  of  the  code  is  an  important  step  of  its  development  because 
numerical  results  must  be  tested  to  ensure  their  accuracy.  Once  the  code  has  been 
validated  it  can  be  used  to  model  thrusters  for  which  no  experiments  have  been  conducted 
yet  and  one  can  have  a  level  of  confidence  in  the  results.  Additionally,  the  code  can  be 
used  to  model  “in-flight”  conditions  that  are  impossible  to  create  in  an  experimental 
environment. 
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6.1 :  Summary  of  Results 


The  results  from  the  validation  and  sensitivity  study  show  that  the  DRACO  code 
can  be  used  to  recreate  a  vacuum  chamber  simulation  using  a  Hall  Thruster.  As  a  result  of 
this  research  the  following  conclusions  can  be  drawn  about  the  use  of  DRACO  to  model 
an  HPHall  source  in  a  vacuum  chamber: 

1.  The  newly  developed  HPHall  source  that  samples  at  the  exit  plane  of  the  thruster 
is  the  best  model  for  a  Hall  Thruster  that  is  currently  available.  The  source 
produces  Faraday  probe  results  that  almost  exactly  match  the  data  from  the 
experiment  conducted  at  AFRL. 

2.  The  Boltzmann  solver  with  the  polytropic  temperature  model  produces  results  that 
are  better  than  the  Boltzmann  solver  that  uses  the  constant  temperature  model  or  a 
simulation  that  does  not  solve  for  the  electric  field. 

3.  Modeling  CEX  collisions  is  important  to  achieve  accurate  simulation  results 
outside  of  the  main  beam  of  the  thruster  for  neutral  densities  on  the  order  of  10  . 

4.  The  collision  model  that  produces  the  best  results  for  this  type  of  simulation  is  the 
MCC  method  using  a  background  neutral  density  that  analytically  predicts  the 
neutral  flow  from  the  thruster  and  projects  that  density  onto  the  background 
neutral  distribution.  This  method  only  models  CEX  collisions  and  does  not  model 
the  actual  neutral  particle  motion  and  as  a  result  is  the  least  computationally 
demanding  of  all  the  collision  methods  studied. 

5.  A  non-uniform  mesh  can  be  implemented  to  produce  results  that  are  as  accurate 
as  a  uniform  mesh,  if  not  more  in  regions  of  finer  meshes.  Additionally,  a  non- 
uniform  mesh  can  be  less  computationally  demanding  because  there  are  fewer 
cells  than  in  a  uniform  mesh. 

6.  A  non-uniform  mesh  can  be  implemented  to  produce  a  more  refined  mesh  at  the 
exit  plane  of  a  thruster  so  that  a  more  accurate  representation  of  the  electric  field 
generated  in  HPHall  can  be  loaded  into  DRACO. 
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6.2:  Suggestions  for  Future  Work 

Although  the  results  of  this  thesis  are  promising  for  the  DRACO  code,  further 
development  and  validation  is  still  required  of  the  code.  The  following  is  a  list  of  future 
work  that  is  still  to  come  for  the  code: 

1.  Further  development  of  the  HPHall  source  to  load  both  a  z-direction  electric  field 
and  also  load  a  radial  electric  field  and  convert  it  into  a  x  and  y  direction  electric 
field.  This  adjustment  to  the  source  could  have  the  potential  to  further  increase  the 
accuracy  of  the  results. 

2.  Run  more  computationally  demanding  simulations  using  parallel  processing  with 
a  larger  domain  and  significantly  more  particles.  The  simulations  in  this  thesis 
contain  approximately  1-3  million  particles  but  a  simulation  run  using  parallel 
processing  could  increase  the  number  of  particles  by  one  or  two  orders  of 
magnitude. 

3.  Investigate  the  use  of  multiple  values  of  specific  heat  ratio  for  use  in  a  Boltzmann 
temperature  model;  possibly  relate  the  specific  heat  ratio  to  plasma  density.  This 
would  be  in  an  attempt  to  improve  the  polytropic  temperature  results  while 
maintaining  the  plasma  potential  results. 

4.  Run  simulations  using  a  better  random  number  generator,  such  as  a  Mersenne 
Twister  RNG,  to  reduce  noise  in  results. 

5.  Run  simulations  using  a  larger  domain  with  the  surface  geometry  of  the  vacuum 
chamber  included.  This  could  be  used  to  maintain  the  neutral  density  for  use  in 
MCC  simulations  that  track  particles  as  well  as  DSMC  simulations. 
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